王貴春，曹宗恒

(鄭州大學(xué) 土木工程學(xué)院 河南 鄭州 450001)

0 引言

渦激振動是由結(jié)構(gòu)尾流區(qū)漩渦脫落引起的結(jié)構(gòu)周期性振動，具有強迫、自激及限幅的特點。斜拉索在流體作用下都可能發(fā)生渦激振動，從而導(dǎo)致結(jié)構(gòu)疲勞破壞，影響結(jié)構(gòu)正常使用。所以研究斜拉索渦激振動響應(yīng)特性具有重要意義。

Feng[1]在1968年進行風洞試驗，研究了單圓柱體渦激振動響應(yīng)。隨著計算機技術(shù)的迅速發(fā)展，很多學(xué)者采用計算流體動力學(xué)(compute fluid dynamics,CFD)方法研究圓柱體渦激振動響應(yīng)特性。李廣望[2]采用任意Lagrange-Euler方法計算了圓柱體的渦激振動響應(yīng)。徐楓[3]采用動網(wǎng)格技術(shù)對不同截面形式的柱體進行了數(shù)值模擬，但是柱體運動引起的計算網(wǎng)格畸變、纏繞等問題不能得到有效解決。劉志文等[4]采用“剛性區(qū)域+動網(wǎng)格區(qū)域+靜止網(wǎng)格區(qū)域”的方式建立網(wǎng)格，對矩形斷面橋梁渦激振動進行計算，計算過程中雖能保證壁面附近計算網(wǎng)格的質(zhì)量，但是運動網(wǎng)格劃分方式較復(fù)雜且不同劃分方式網(wǎng)格間的連接容易出現(xiàn)問題，尤其是當網(wǎng)格出現(xiàn)變形時，動網(wǎng)格區(qū)域計算網(wǎng)格質(zhì)量大幅下降，容易影響計算結(jié)果的精度。劉軍[5]采用Fluent對高雷諾數(shù)下輸電導(dǎo)線渦激振動進行了分析，并詳細分析了高質(zhì)量比和低質(zhì)量比參數(shù)對渦激振動的影響，其研究結(jié)果表明：低質(zhì)量比渦激振動振幅要大于高質(zhì)量比渦激振動振幅。曲寶旭[6]采用重疊網(wǎng)格法對具有彈性支撐不同自由度下輸電導(dǎo)線的渦激振動進行了分析，結(jié)果表明：不同自由度對渦激振動橫風向振動具有一定影響。

當前對圓柱渦激振動數(shù)值模擬主要是采用動網(wǎng)格技術(shù)，采用重疊網(wǎng)格法的相關(guān)報道相對較少。本文采用重疊網(wǎng)格技術(shù)建立網(wǎng)格，應(yīng)用Newmark-β法求解結(jié)構(gòu)振動方程，并且編制UDF(user defined functions)，嵌入到Fluent中，通過宏命令DEFINE_ZONE_MOTION傳遞運動速度，以更新運動網(wǎng)格位置，對斜拉索結(jié)構(gòu)渦激振動響應(yīng)進行數(shù)值模擬，從而解決結(jié)構(gòu)振動過程中可能產(chǎn)生的網(wǎng)格畸變和負體積問題，實現(xiàn)流固耦合計算。

1 理論基礎(chǔ)

1.1 流體控制方程

不可壓縮流體N-S方程和連續(xù)性方程分別為

(1)

(2)

式中：υ=μ/ρ為分子黏度，μ為動力黏度；fi為單位質(zhì)量流體所受到的體積力(質(zhì)量力)；p為流體介質(zhì)微元體上的壓力；ρ表示流體密度；ui表示x、y、z方向的流體速度分量。

1.2 結(jié)構(gòu)振動控制方程

以文獻[7]中斜拉索實驗?zāi)Ｐ蜑楸尘斑M行斜拉索渦激振動分析，斜拉索直徑D設(shè)置為0.066 m。如圖1所示，斜拉索可簡化為兩自由度質(zhì)量-彈簧-阻尼系統(tǒng)，振動方程表示為

圖1 兩自由度質(zhì)量-彈簧-阻尼體系

(3)

2 數(shù)值模型建立

在劃分網(wǎng)格之前，要確定計算域的大小，合理的計算域大小可以使數(shù)值計算效率提高，也能保持計算結(jié)果的準確性。通?？梢愿鶕?jù)阻塞率的大小來確定計算域大小，在CFD模型中阻塞率一般控制在5%以內(nèi)。二維流場計算域大小為24D×50D，其中D為圓柱的直徑，如圖2所示。定義圓柱截面形心所在位置為計算域坐標原點，設(shè)定重疊區(qū)域直徑大小為3D，坐標原點到入口邊界(INLET)的水平距離為15D，與出口邊界(OUTLET)的水平距離為30D。上下兩側(cè)均為對稱邊界(SYMMETRY)，其到坐標原點的垂直距離均為12D。拉索表面為無滑移壁面邊界(NO-SLIPWALL)。在進行拉索渦激振動計算時，湍流模采用剪切應(yīng)力輸運模型(shear stress transition，SST)k-ω模型，湍流度設(shè)置為5%。該模型能較好地處理壁面附近的流動得到可靠的結(jié)果，適用于存在逆壓梯度和流動分離的圓柱繞流場計算[8]。

圖2 重疊網(wǎng)格劃分

2.1 重疊網(wǎng)格劃分及邊界條件

重疊網(wǎng)格由背景網(wǎng)格和前景網(wǎng)格組成，也即流體計算域根據(jù)計算需要被劃分為多個部件區(qū)域。首先在各區(qū)域單獨創(chuàng)建網(wǎng)格，然后通過插值計算實現(xiàn)相鄰計算域間的流場數(shù)據(jù)交換。與傳統(tǒng)動網(wǎng)格計算相比，重疊網(wǎng)格在計算過程中，不會產(chǎn)生網(wǎng)格纏繞、畸變和負體積等問題，能始終保持計算網(wǎng)格的質(zhì)量，從而保證計算結(jié)果的準確性。

采用CFD軟件Fluent中的重疊網(wǎng)格(overset mesh)技術(shù)，進行圓柱體兩自由度渦激振動響應(yīng)計算。重疊網(wǎng)格技術(shù)的實施，包括兩個主要步驟：1)在不同網(wǎng)格間建立通信關(guān)系；2)利用從步驟1中獲得的信息，進行數(shù)值求解。

2.2 網(wǎng)格獨立性與時間無關(guān)性檢驗

表1 網(wǎng)格方案詳細參數(shù)設(shè)置

從圖3(a)可以看出在不同網(wǎng)格方案設(shè)置下，拉索氣動響應(yīng)隨著時間步的減小都能達到一個穩(wěn)定的值，不同的是網(wǎng)格數(shù)量越大以及寬高比越小，也就是網(wǎng)格的整體質(zhì)量越好(A1,A2網(wǎng)格)，得到的阻力系數(shù)均值較A3網(wǎng)格方案得到的結(jié)果稍大。從圖3(b)中可以看出，在不同網(wǎng)格方案下，采用A1和A2網(wǎng)格得到的St數(shù)在整體趨勢上更加吻合，結(jié)果也比較接近?？梢缘贸鼋Y(jié)論，采用A1和A2網(wǎng)格方案，得到St數(shù)與其他學(xué)者的結(jié)果比較一致，如表2所示。可以認為A2、A3網(wǎng)格方案已通過網(wǎng)格獨立性以及時間無關(guān)性檢驗，繼續(xù)對網(wǎng)格進行加密意義不大。綜上所述，在考慮計算精確度和計算效率下，選擇A2網(wǎng)格方案和時間步大小為0.002 s進行數(shù)值模擬。

表2 Re=3 900時St數(shù)對比

圖3 阻力系數(shù)均值與St數(shù)隨時間步大小的變化

3 斜拉索渦激振動響應(yīng)特性分析

拉索在實際服役時，渦激振動起振的折減風速為Ur=U/(fn·D)，通常小于10，所以選取的計算折減風速Ur范圍為3.5～7.5，相應(yīng)的雷諾數(shù)Re位于6 500～14 000之間。

3.1 St數(shù)及渦脫頻率分析

圖4(a)中可以看出，在折減風速Ur=3.5～4.64之間時，St數(shù)隨折減風速變化呈線性下降的趨勢；當Ur=4.727時，St=0.203；在Ur=4.64～4.94之間時，隨著Ur的增大St呈倒“V”字形變化，即St數(shù)先突然增大，又急劇減小，此時對應(yīng)拉索發(fā)生渦激振動的折減風速區(qū)間，可以推斷這個區(qū)間內(nèi)流場變化比較劇烈；當Ur>4.94時，St數(shù)又突然增大之后呈線性減小趨勢。如圖4(b)所示，隨著折減風速增大，由于拉索發(fā)生了頻率鎖定現(xiàn)象(lock-in)，所以在某一風速區(qū)間內(nèi)拉索渦脫頻率并沒有線性增大；當渦脫頻率約為3.494 Hz時，對應(yīng)的折減風速為4.94；當Ur>4.59時，拉索渦脫頻率沒有隨著折減風速的增大而增大，而是與在靠近拉索的一階固有頻率，所以，可以認為此處為鎖定區(qū)間的起點；當Ur=5.04時，拉索渦脫頻率與結(jié)構(gòu)固有頻率解鎖，渦脫頻率隨折減風速的增加開始呈線性變化。從圖4(b)中可以看出，拉索發(fā)生鎖定的區(qū)間比較窄，也就說，拉索渦激振動對風速的變化非常敏感。

圖4 不同折減風速下渦激振動St數(shù)與升力系數(shù)頻率分析結(jié)果

3.2 振幅響應(yīng)分析

圖5給出了不同折減風速下拉索橫風向無量綱位移時程曲線。從圖5中可以看出，在高質(zhì)量比下拉索在折減風速較低時，位移呈現(xiàn)微幅簡諧振動，拉索的振動幅度也很小。隨著折減風速的增加拉索振動幅度逐漸增大，并且出現(xiàn)了“拍振”現(xiàn)象。橫向位移在折減風速Ur=4.67～4.94之間顯著增加。

圖5 典型折減風速下拉索橫向振動位移時程曲線

如圖6(a)所示，可以看出，在鎖定區(qū)間拉索橫風向振動頻率接近拉索固有頻率。在折減風速Ur=4.64～4.94之間時，漩渦脫落頻率與拉索固有頻率值非常接近，即拉索發(fā)生了“鎖定”，拉索作大幅振動，并且在折減風速Ur=4.94時，拉索達到最大振幅0.20D左右(如圖6(b)所示)。

圖6 不同折減風速下拉索橫向主振頻率與最大位移

3.3 拉索渦激振動軌跡分析

圖7(a)～(g)給出了不同折減風速下拉索渦激振動的軌跡圖。拉索渦激振動是以空氣為介質(zhì)進行計算，拉索的質(zhì)量比m*=159.4，此時拉索應(yīng)屬于高質(zhì)量比下渦激振動。高質(zhì)量比情況下，渦激振動的運動軌跡形式更加復(fù)雜。在較低或者較高折減風速下，拉索振動軌跡主要呈橢圓形式，同時運動的方向和位置也在不斷發(fā)生變化。在接近結(jié)構(gòu)固有頻率時，拉索在鎖定區(qū)內(nèi)也即流固耦合作用較強時，拉索運動軌跡大致呈“8”字形。低質(zhì)量比圓柱渦激振動軌跡一般呈單一的“8”字形，高質(zhì)量比結(jié)構(gòu)渦激振動的運動形式更加豐富，在高質(zhì)量比參數(shù)下拉索不再是單一的“8”字形運動，這與文獻[6]的結(jié)論一致。

圖7 拉索渦激振動質(zhì)心運動軌跡

4 質(zhì)量比對渦激振動影響分析

本節(jié)設(shè)定質(zhì)量比m*=7.9作為低質(zhì)量比情況，阻尼比保持ζ=0.546%不變(相對應(yīng)的質(zhì)量阻尼m*ζ=0.041 5)，對拉索進行數(shù)值模擬。

從圖8中可以看到，在質(zhì)量比為159.4時，升阻力系數(shù)的均值是比較接近的，也就是說，在高質(zhì)量比參數(shù)下，流場的變化并不是很強烈。對比不同質(zhì)量比下拉索升阻力系數(shù)均值可以發(fā)現(xiàn)，低質(zhì)量比參數(shù)下拉索升阻力系數(shù)均值變化的范圍更廣，橫風向振動大小對升阻力大小有較大的影響，在結(jié)構(gòu)振幅較大時，其升阻力也有一個非常顯著的增大過程，但隨著振幅的降低，其升阻力又重新降低，直到保持一個穩(wěn)定的值?？梢酝茢?，低質(zhì)量比參數(shù)下流場與結(jié)構(gòu)之間的相互作用也更加顯著。

圖8 不同折減風速及質(zhì)量比下升阻力系數(shù)均值對比

如圖9、10所示，相較于高質(zhì)量比情況，拉索渦激振動的鎖定區(qū)間也有所變大，由鎖定區(qū)過渡到非鎖定區(qū)較慢。當頻率比為1時，拉索達到最大振幅為0.82D，為高質(zhì)量比參數(shù)下最大振幅的4.1倍。在設(shè)置計算的有限工況內(nèi)，渦激振動處于鎖定區(qū)間時的鎖定頻率比在1.02左右，也就是說，鎖定頻率發(fā)生了小幅度的偏移。當頻率比達2.318時，渦激振動振幅顯著減小至0.05D，此時渦激振動已脫離鎖定區(qū)間。

圖9 不同折減風速下橫風向最大位移變化

對比圖8～10可以推斷：質(zhì)量比越小，結(jié)構(gòu)發(fā)生渦激振動時鎖定區(qū)間越寬，最大振幅也越大。隨著質(zhì)量比減小，發(fā)生鎖定時的頻率比也會有一定的偏移。

圖10 不同折減風速下頻率比變化

5 結(jié)論

通過對斜拉索進行渦激振動數(shù)值模擬，得出了4個結(jié)論。1)拉索發(fā)生渦激振動時，阻力系數(shù)均值首先降低至最小值，而后又顯著增加到最大值。之后隨著折減風速的繼續(xù)增大，阻力系數(shù)均值呈不斷減小趨勢。升力系數(shù)均方根值隨折減速度的增大(Ur=3.5～4.51)基本保持不變，之后升力系數(shù)突然降低，又在Ur=4.59～4.94之間保持顯著增大的趨勢，直至最大值，之后隨折減風速的增加呈減小趨勢；2)對于高質(zhì)量比渦激振動，其振幅雖然較低，但由于振動產(chǎn)生的疲勞破壞要特別注意。只有在接近拉索一階固有頻率的范圍內(nèi)，拉索振幅顯著增大，達到了0.2D，之后雖折減風速進一步增加，但振幅迅速減小，拉索振幅的變化對于折減風速的變化非常敏感；3)對拉索渦激振動的運動軌跡進行了詳細的分析，拉索屬于高質(zhì)量比渦激振動，其運動軌跡更加豐富，只有在接近發(fā)生渦激振動的范圍內(nèi)，拉索大致呈“8”字形運動軌跡，這與低質(zhì)量比圓柱的渦激振動軌跡有顯著差異；4)質(zhì)量比會對拉索渦激振動產(chǎn)生一定的影響，是影響渦激振動的重要因素。通過對低質(zhì)量比渦激振動分析，可以得出結(jié)論：質(zhì)量比越小，振動鎖定區(qū)間越大，拉索在鎖定區(qū)間內(nèi)振幅也越大。