有攻角方柱渦致振動的數值研究

湯兆烈， 周本謀

(南京理工大學 瞬態物理重點實驗室，南京 210094)

流致振動廣泛存在于各類自然環境與工程領域中，例如系泊電纜，海洋鉆井立管，輸油立管等。振動會使結構疲勞發生破壞，但是隨著研究進展，人們發現也可以利用這種現象獲取能源[1-4]。結構流致振動的形式比較多，渦致振動和馳振是工程中常見的兩個現象。對于通常形狀的結構物而言，渦激振動是必然伴隨的現象，而馳振則因鈍體截面的形狀不同而有所差異，多發生于具有非流線型或非對稱截面的結構。其中圓柱作為渦致振動典型例子得到廣泛研究[5]，同時方柱[6]作為馳振代表，矩形柱[7]、三角柱[8]往往會是二者的結合。

與無攻角方柱或矩形柱相比，有攻角方柱研究較少。圓柱渦致振動的響應可以分成三個區間[9]：初始分支(initial branch)、上部分支(upper branch)和下部分支(lower branch)。Nemes等[10]實驗研究了攻角α對方柱流致振動的影響，發現當攻角在0°～10°馳振占主導；當攻角在25°以上時，渦激振動占主導；當時攻角在10°～25°時，振動響應變得復雜，渦激振動與馳振共同起作用，出現了一個“高分支”(Higher Branch,HB)，該分支的振幅比上支還要高。后續Zhao等[11]選取攻角α=20°，更加精確和系統的描述這種現象，發現高分支是一種高次諧波，是由于對稱性被破壞導致的，同時也能解釋對稱的圓形截面不會出現這種現象。其他人不同截面也同樣得到了高次諧波，如三角形截面[12]、橢圓截面[13]。數值實驗方面，Zhao等[14-15]都在攻角22.5°得到了Leontini等[16]使用譜元法計算攻角45°時低雷諾數下不同程度圓角對方柱渦致振動的影響。其中質量比固定為2，當沒有圓角或者圓角比較小時，隨著折合速度或系統剛度的變化，他們也得到了高次諧波現象的出現和消失。這預示著HB的出現影響因素很多。不過上述試驗和數值模擬中，僅圍繞同一種質量，并未考慮質量比的影響。而根據現有的渦激振動研究可知，系統剛度、質量比、雷諾數、阻尼都是是影響渦致振動的關鍵性因素。Govardhan等[17]發現圓柱低質量比鎖定區間比高質量比時寬。Sen等[18]研究了質量比對無攻角方柱渦致振動的影響，發現隨著質量比增大，馳振開始出現的雷諾數下降。練繼建等[19]實驗研究了質量比和系統剛度對不同攻角下方柱渦激振動的影響，他們僅僅是關心振幅和頻率，并沒有提到次諧波現象。為此，本文利用譜元法研究不同質量比以及系統剛度對有攻角方柱振動的影響。本文的研究目的：研究不同質量比以及系統剛度對方柱渦致振動的影響，尤其是高次諧波出現的質量比以及系統剛度的范圍。

1 數值方法及驗證

1.1 控制方程

方柱與來流形成一定的攻角α，并且隨著渦的脫落發生振動，由于發生渦致振動的柱體在流向上的振幅遠小于橫向振幅，因此本文僅對橫向振動進行研究，模型如圖1所示。參照Zhao等的研究，選取攻角α=20°來進行計算。

圖1 方柱渦致振動模型Fig.1 Sketch of VIV of a square cylinder

該問題的控制方程是無量綱二維不可壓黏性N-S方程以及物體運動方程。

無量綱N-S方程

(1)

(2)

無量綱物體運動方程

(3)

1.2 數值方法

計算區域為：入口邊界以及上下側邊界距方柱中心為20D，出口邊界距方柱中心為30D。入口邊界以及上下側邊界設定為u=1,v=0，出口邊界設定為壓力出口，柱體表面設為無滑移壁面條件。方程的求解采用了譜元法(Spectral-Element Method)，使用的求解軟件是Nektar++[20]。空間區域分成了457個四邊形單元，如圖2(a)所示，形函數采用8階Lagrange多項式，插值點采用Gauss-Lobatto-Legendre點，計算點如圖2(b)所示。初始條件都設為零，根據Navrose等[21]，本文計算的質量比，不同的初始條件對結果幾乎沒影響。

圖2 計算網格Fig.2 Computational mesh

譜元法的收斂性驗證是通過保持大的網格不變，通過增加多項式形函數的階數以達到收斂，這樣的過程稱為p收斂。選取m*=2,Ur=3，通過增加網格內插值點的數量達到收斂，結果如表1所示。由表可以看出，p=8時，與更高階數的差別在0.5%以內，因此，后續計算均選用p=8。

表1 插值階數收斂性驗證Tab.1 Convergence with different polynomial orders

1.3 結果驗證

與Zhao等的研究對比，發現除了5≤Ur≤7以及20，本文計算的其它折合速度下橫向物體位移幅Ay差別在5%以內，Zhao等研究中5≤Ur≤7柱體位移隨時間的變化是周期性變化的，而本文計算結果不是周期性的。 Leontini等研究中計算條件為α=45°質量比m*=3，Re=200， 他們的計算結果顯示3.5≤Ur≤7.8柱體位移隨時間的變化都是非周期性的。另外Nemes等的實驗結果也表明，上部分支中柱體位移隨時間變化是非周期性。而折合速度Ur=20時，柱體位移振幅變化太劇烈，導致和文獻對比存在一定的差別。

圖3 與Zhao等研究中數值計算結果的對比Fig.3 Comparison between present results with the numerical result by Zhao et al

2 結果及分析

2.1 質量比的影響

首先研究質量比的影響，計算條件為Re=200選， 選取三個不同質量比m*=2,m*=5,m*=10， 折合速度范圍Ur=2～20， 振幅響應如圖4所示。

圖4 不同質量比下柱體橫向振幅Fig.4 Amplitude in the cross-flow direction with different mass ratios

圖5是不同質量比下柱體位移以及升力系數功率圖譜，可以從中看到隨著折合速度的變化，存在一個階段位移和升力系數時程曲線是高次諧波。

圖5中，左側是物體橫向位移功率譜云圖，右側是升力系數功率譜云圖。具體的做法是在某一給定折合速度下將物體橫向位移或升力系數隨時間的變化曲線利用快速傅里葉變換(Fast Fourier Transformation，FFT)得到功率譜(頻率除以U/D無量綱化)，將結果都除以最大功率歸一化，然后將不同折合速度得到的結果整合在一起形成云圖。云圖顏色的深淺代表著某一折合速度下，某一頻率的功率的均一化之后的值。

結合圖4和圖5，可以看到與圓柱渦致振動類似，隨著折合速度的增加，有攻角方柱振動響應也可以分

圖5 不同質量比下橫向位移、升力系數功率譜Fig.5 Power spectrum contour plots of body oscillation, lift coefficient with different mass ratios

為起始支、上部分支以及下部分支。不同的是，低質量比m*=2,m*=5時，出現了Nemes等描述的高分支，該分支的特點是振幅和升力系數隨時間變化是周期性非正弦，且進行傅里葉變換得到的基波頻率整數倍的各次分量。隨著質量比變為m*=10，這一分支消失。可以看到有攻角方柱渦致振動響應不僅跟折合速度有關，還跟質量比有很大的關系。為了進一步說明這一點，固定折合速度Ur=10， 質量比從2變化到10，結果如圖6所示。當質量比大于6時，高分支消失。

圖6 Ur=10不同質量比物體位移、升力系數功率圖Fig.6 Power spectrum contour plots at Ur=10

2.2 有效剛度

為了更加清楚地描述高分支，以m*=2,Ur=9為例。圖7是相應柱體橫向位移和升力系數的時程曲線以及功率圖譜。柱體靜止時St=0.154，從功率圖譜中可以看到位移和升力系數的頻率主要集中在St, 1/2St附近。

圖7 m*=2, Ur=9物體位移升力 系數歷史曲線及功率譜Fig.7 Time traces and normalized power spectra of displacement, lift coefficient at m*=2, Ur=9

圖8是一個周期內不同時刻柱體尾跡渦量圖。從渦量圖中可看出，一個周期內交替脫落兩組渦，Nemes等稱之為2(2S)模式。可以看到，渦的強度是不同的，類似于P+S。造成這種渦脫落形式的主要原因是流體剪切層和柱體后緣之間相互作用。柱體尾部在一定情況下由于運動會形成后緣渦，同時在柱體的另一邊會形成一個與之相反的前緣渦，這兩個渦交替脫落共同造成2(2S)模式。

圖8 m*=2, Ur=9渦量圖Fig.8 Vorticity contour at m*=2, Ur=9

為了描述圓柱渦致振動“鎖定”階段，Shiels等[22]提出了有效剛度概念，具體的做法是采用另一種無量綱運動方程

(4)

式中：b*=b/(ρUD);k*=k/(ρU2)。

在本文計算中不考慮阻尼的影響，因此上式變為

(5)

然后假定鎖定階段物體橫向位移和升力系數都是正弦變化，并且他們的相位差為0或π，此時，位移幅值以及升力系數幅值滿足一定的關系，它們之間的比值稱為有效剛度。利用有效剛度的思想，下面同樣假定高次諧波出現的范圍內，物體橫向位移以及升力系數滿足下面的形式

Y=A0+A1sin(ωt+φ1)+A2sin(2ωt+φ2)+ …+Ansin(nωt+φn)

(6a)

Cl=C0+C1sin(ωt+φ1)+C2sin(2ωt+φ2)+ …+Cnsin(nωt+φn)

(6b)

式中：ω=2πf，f為基本頻率。

將物體橫向位移以及升力系數表達式代入式(5)中得到

k*A0+A1(k*-m*ω2)sin(ωt+φ1)+A2(k*-m*ω2)sin(2ωt+φ2)+…+An(k*-m*ω2)sin(nωt+φn)= 1/2[C0+F1sin(ωt+φ1)+C2sin(2ωt+φ2)+…+Cnsin(nωt+φn)]

(7)

為了使式子成立，則

A0k*=1/2C0A1(k*-m*ω2)=1/2C1A2(k*-4m*ω2)=1/2C2……An(k*-n2m*ω2)=1/2Cn

(8)

在出現高次次諧波這一階段，可以認為系統是線性的。另外，沒有兩組(k*,m*)能夠滿足所有的k*-n2m*ω2相同，這說明系統剛度和質量比這兩個量在描述HB這一階段是獨立的。因此，采用折合速度作為參量是不能描述清楚的。

現令

(9)

稱它們為第i有效剛度，隨著有效剛度次數增大，相應的值也越來越大，這意味著升力系數與物體振幅之間的比值也越來越大，這也能解釋為什么物體橫向位移功率譜中高頻率成分難以被發現，而升力系數高頻率成分依然很明顯。

2.3 k*-m*圖像

由“2.2”節得知，高分支有多個有效剛度，為了研究清楚這些有效剛度對結果的影響，選取了一系列系統剛度和質量比作為計算點，結果如圖9所示。

圖9中，正方形表示高分支，三角形表示上部分支，圓形表示下部分支。高分支出現的范圍大致是

k*-4π2×0.142m*≤0k*-π2×0.142m*≤2k*-π2×0.142m*≥0.2k*-4π2×0.142m*≥-3.2

(10)

k-π2×0.142m*≥0.2，k-π2×0.142m2≥0.2對應著

圖9 k*-m*計算結果Fig.9 Results with different k* and m*

Ur=7,14， 這和“2.1”節得到的高分支出現最大Ur區間相吻合。 另外，當m*≥7，沒有發現高分支，稱該質量比為臨界質量比。

發現了一些不同形式的高分支的質量比和剛度(m*,k*)點， 雖然位移時程曲線是高次諧波，但基本頻率是1/3St, 1/4St， 如圖10和圖11所示。

圖10 m*=6, k*=4物體位移時程曲線及功率譜Fig.10 Time traces and normalized power spectra of displacement at m*=6, k*=4

圖11 m*=4.5, k*=2.5物體位移時程曲線及功率譜Fig.11 Time traces and normalized power spectra of displacement at m*=4.5, k*=2.5

更低基本頻率的發現，更加說明了有攻角方柱渦致振動的復雜性。

和圓柱渦致振動一樣，“拍”現象也被發現(見圖12)。“拍”是渦致振動中一種不太常見的物理現象，產生的必要條件是兩個振動頻率比較接近，這時候位移響應和升力系數的值表現出調諧狀態。

圖12 m*=2, k*=0.5出現的“拍”現象Fig.12 Beating at m*=2, k*=0.5

不足之處，沒有考慮雷諾數、攻角α以及系統阻尼的影響。可預見的是雷諾數和攻角α影響的是渦脫落頻率，進而影響有效剛度。而系統阻尼影響的是物體運動與升力之間的相位差，以及最大振幅與升力系數之間的比值。

3 結 論

本文使用譜元法計算了有攻角方柱橫向自由振動，重點研究了質量比以及系統剛度對振動響應特征的影響。現得到的結論如下：

隨著質量比的增大，高分支出現的區間逐漸變小直到消失。比起折合速度，有效剛度更能描述這一階段。與圓柱不同的是，圓柱一個有效剛度就可以，而且當有效剛度相等時，他們的響應相同；而對于與來流呈一定角度的方柱而言，高分支有多個有效剛度，因此不會有兩組(k*,m*)，使得振動響應完全相同。高分支出現的范圍主要與第一、第二剛度有關。另外，發現高分支的基本頻率不僅僅是1/2St，還可以是1/3St甚至1/4St。