有攻角方柱渦致振動(dòng)的數(shù)值研究

湯兆烈， 周本謀

(南京理工大學(xué) 瞬態(tài)物理重點(diǎn)實(shí)驗(yàn)室，南京 210094)

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

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

1 數(shù)值方法及驗(yàn)證

1.1 控制方程

方柱與來(lái)流形成一定的攻角α，并且隨著渦的脫落發(fā)生振動(dòng)，由于發(fā)生渦致振動(dòng)的柱體在流向上的振幅遠(yuǎn)小于橫向振幅，因此本文僅對(duì)橫向振動(dòng)進(jìn)行研究，模型如圖1所示。參照Z(yǔ)hao等的研究，選取攻角α=20°來(lái)進(jìn)行計(jì)算。

圖1 方柱渦致振動(dòng)模型Fig.1 Sketch of VIV of a square cylinder

該問(wèn)題的控制方程是無(wú)量綱二維不可壓黏性N-S方程以及物體運(yùn)動(dòng)方程。

無(wú)量綱N-S方程

(1)

(2)

無(wú)量綱物體運(yùn)動(dòng)方程

(3)

1.2 數(shù)值方法

計(jì)算區(qū)域?yàn)椋喝肟谶吔缫约吧舷聜?cè)邊界距方柱中心為20D，出口邊界距方柱中心為30D。入口邊界以及上下側(cè)邊界設(shè)定為u=1,v=0，出口邊界設(shè)定為壓力出口，柱體表面設(shè)為無(wú)滑移壁面條件。方程的求解采用了譜元法(Spectral-Element Method)，使用的求解軟件是Nektar++[20]?？臻g區(qū)域分成了457個(gè)四邊形單元，如圖2(a)所示，形函數(shù)采用8階Lagrange多項(xiàng)式，插值點(diǎn)采用Gauss-Lobatto-Legendre點(diǎn)，計(jì)算點(diǎn)如圖2(b)所示。初始條件都設(shè)為零，根據(jù)Navrose等[21]，本文計(jì)算的質(zhì)量比，不同的初始條件對(duì)結(jié)果幾乎沒(méi)影響。

圖2 計(jì)算網(wǎng)格Fig.2 Computational mesh

譜元法的收斂性驗(yàn)證是通過(guò)保持大的網(wǎng)格不變，通過(guò)增加多項(xiàng)式形函數(shù)的階數(shù)以達(dá)到收斂，這樣的過(guò)程稱為p收斂。選取m*=2,Ur=3，通過(guò)增加網(wǎng)格內(nèi)插值點(diǎn)的數(shù)量達(dá)到收斂，結(jié)果如表1所示。由表可以看出，p=8時(shí)，與更高階數(shù)的差別在0.5%以內(nèi)，因此，后續(xù)計(jì)算均選用p=8。

表1 插值階數(shù)收斂性驗(yàn)證Tab.1 Convergence with different polynomial orders

1.3 結(jié)果驗(yàn)證

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

圖3 與Zhao等研究中數(shù)值計(jì)算結(jié)果的對(duì)比Fig.3 Comparison between present results with the numerical result by Zhao et al

2 結(jié)果及分析

2.1 質(zhì)量比的影響

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

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

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

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

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

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

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

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

2.2 有效剛度

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

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

圖8是一個(gè)周期內(nèi)不同時(shí)刻柱體尾跡渦量圖。從渦量圖中可看出，一個(gè)周期內(nèi)交替脫落兩組渦，Nemes等稱之為2(2S)模式?？梢钥吹?，渦的強(qiáng)度是不同的，類似于P+S。造成這種渦脫落形式的主要原因是流體剪切層和柱體后緣之間相互作用。柱體尾部在一定情況下由于運(yùn)動(dòng)會(huì)形成后緣渦，同時(shí)在柱體的另一邊會(huì)形成一個(gè)與之相反的前緣渦，這兩個(gè)渦交替脫落共同造成2(2S)模式。

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

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

(4)

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

在本文計(jì)算中不考慮阻尼的影響，因此上式變?yōu)?/p>

(5)

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

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為基本頻率。

將物體橫向位移以及升力系數(shù)表達(dá)式代入式(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)

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

現(xiàn)令

(9)

稱它們?yōu)榈趇有效剛度，隨著有效剛度次數(shù)增大，相應(yīng)的值也越來(lái)越大，這意味著升力系數(shù)與物體振幅之間的比值也越來(lái)越大，這也能解釋為什么物體橫向位移功率譜中高頻率成分難以被發(fā)現(xiàn)，而升力系數(shù)高頻率成分依然很明顯。

2.3 k*-m*圖像

由“2.2”節(jié)得知，高分支有多個(gè)有效剛度，為了研究清楚這些有效剛度對(duì)結(jié)果的影響，選取了一系列系統(tǒng)剛度和質(zhì)量比作為計(jì)算點(diǎn)，結(jié)果如圖9所示。

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

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對(duì)應(yīng)著

圖9 k*-m*計(jì)算結(jié)果Fig.9 Results with different k* and m*

Ur=7,14， 這和“2.1”節(jié)得到的高分支出現(xiàn)最大Ur區(qū)間相吻合。 另外，當(dāng)m*≥7，沒(méi)有發(fā)現(xiàn)高分支，稱該質(zhì)量比為臨界質(zhì)量比。

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

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

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

更低基本頻率的發(fā)現(xiàn)，更加說(shuō)明了有攻角方柱渦致振動(dòng)的復(fù)雜性。

和圓柱渦致振動(dòng)一樣，“拍”現(xiàn)象也被發(fā)現(xiàn)(見(jiàn)圖12)。“拍”是渦致振動(dòng)中一種不太常見(jiàn)的物理現(xiàn)象，產(chǎn)生的必要條件是兩個(gè)振動(dòng)頻率比較接近，這時(shí)候位移響應(yīng)和升力系數(shù)的值表現(xiàn)出調(diào)諧狀態(tài)。

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

不足之處，沒(méi)有考慮雷諾數(shù)、攻角α以及系統(tǒng)阻尼的影響?？深A(yù)見(jiàn)的是雷諾數(shù)和攻角α影響的是渦脫落頻率，進(jìn)而影響有效剛度。而系統(tǒng)阻尼影響的是物體運(yùn)動(dòng)與升力之間的相位差，以及最大振幅與升力系數(shù)之間的比值。

3 結(jié) 論

本文使用譜元法計(jì)算了有攻角方柱橫向自由振動(dòng)，重點(diǎn)研究了質(zhì)量比以及系統(tǒng)剛度對(duì)振動(dòng)響應(yīng)特征的影響?，F(xiàn)得到的結(jié)論如下：

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