豎立圓柱渦激振動尖端效應的數值模擬研究

丁 鈺，于哲峰，毛玉明，王吉飛，舒忠平

(1.上海交通大學 航空航天學院，上海 200240;2.上海宇航系統工程研究所，上海 201109)

流體流過柱狀結構的過程是一種非定常的分離流動，位于柱狀體兩側泄落的旋渦會使得結構同時受到橫流向和順流向的脈動壓力，該脈動壓力會引起結構在橫流向和順流向的振動。同時，柱狀體的振動反過來又會改變流場的流動形態以及尾流結構[1]。這種流體-結構相互作用的過程被稱作“渦激振動(vortex-induced vibration，VIV)”。由于結構響應和流體流動特性的高度耦合使得渦激振動成為一種極其復雜的流固耦合過程。其中涉及許多科學上的難題，如湍流、流動分離、剪切層的不完全轉捩以及分離點的漂移等，至今學術界尚無法完全理解和徹底解決這一難題[2]。在實際工程應用中，渦激振動引起的振幅劇烈增加可能會導致材料的疲勞甚至結構的破壞。因此，對于渦激振動問題的研究有著十分重要的學術意義和工程應用意義。

目前對于渦激振動問題的研究，根據流動介質的不同，主要分為兩個領域。其一是對于深海立管等圓柱結構受海水作用下的水動力特性問題[3-4]。另一領域是橋梁建筑物等柱狀體受風載荷作用的響應問題[6-7]。在水和空氣這兩種不同的流動介質作用下，渦激振動的發生都可能會對結構產生嚴重的破壞或者疲勞損傷。關于渦激振動響應的影響因素有很多，比如長徑比、流固質量比、阻尼比、雷諾數以及結構的支承方式等[8-13]都會對流場的流動特性以及結構的振動響應產生很大的影響。同時還有一些學者研究了不同結構外形下流場的流動特性，比如錐形圓柱、波浪形圓柱等[14-15]。

大多數學者使用試驗方法或者數值模擬的方式研究了固定圓柱或者彈性支承的無限長圓柱的繞流流動情況，對于有限長懸臂圓柱的渦激振動自由端效應研究較少。而在實際工程中，往往很多結構都是類似于有限長柱體，例如高聳建筑物、燈柱等。有限長圓柱在自由端處的流動下洗特性對流動的三維特性會產生很大的影響，這和無限長圓柱的繞流流動有很大的區別。最早的有限長圓柱繞流流場模型是Taneda[16]在Re=39～75的條件下通過尾流可視化試驗得到的，該模型提出了復雜的旋渦結構，在試驗里可以清楚地看見流場的下洗運動及其作用效果。但是該模型沒有解釋圓柱附近漩渦動力的詳細情況，也沒有解釋固定圓柱自由端附近的流場。雖然Etzold等[17]也在自己的試驗中得到了在自由端前緣產生的頂渦，并詳細介紹了自由端的流動情況。但是在他們的試驗中所使用的模型都是固定圓柱，所以他們并沒有考慮結構振動對流場特性的影響。后來很多學者針對不同的雷諾數和長徑比對流場的流動特性的影響進行了探索，這些流動模型雖然得到了詳細的漩渦結構及流場特性。但是也都采用的是固定圓柱，忽略了結構產生的位移對流場的影響。另外，對于風載荷作用下的渦激振動的研究，人們大多只考慮了在流體的渦泄頻率和圓柱的固有頻率相接近時的渦激振動現象，針對相對較高風速下的渦激振動響應特性研究較少。雖然一些學者通過針對圓柱模型的試驗研究發現，當風速高于普通渦激振動風速時，渦激振動也可能發生，但其發生的機理尚未明確。比如對于二維圓柱，Durgin等[18]做了一個風洞實驗來研究圓柱的橫風響應。他們得到了普通渦激振動的響應，同時當風速提高至普通渦激振動風速的3倍時，也得到了一種渦激振動響應。根據他們的分析，這種振動的機理被解釋為由于卡門渦街脫落是引起的次頻振動。關于引起結構次頻振動的具體原因也是不得而知。另外，他們都是僅研究了二維流動模型，忽略了流場的三維流動特性。在Matsumoto等[19-20]研究中考慮三維流動特性，他們分別使用均勻的圓柱形拉索和圓臺狀結構模型進行試驗，測量了模型后沿展向不同位置的風速波動。他們的研究結果表明通過對數據進行頻譜分析，發現卡門渦街脫落頻率沿展向變化。則他們認為在雙自由度三維模型中的三維旋渦泄落是導致當風速大于普通渦激振動風速時也產生渦激振動的原因，但是并沒有給出圓柱后方的渦泄頻率沿展向變化的原因。同時，Kawai[21]利用錐形圓柱體模型進行了風洞實驗，得到了普通渦激振動的響應，以及風速為普通渦激振動風速的2.5倍時的響應峰值。在橫風響應的功率譜中，出現了與卡門渦街脫落不同的氣動力峰值。他推斷，當風速在普通渦激振動的2.5倍的情況下，這種氣動力會導致響應峰值。但是沒有確切的證據表明這一點。

近20年來，隨著計算流體動力學(computational fluid dynamics，CFD)理論以及計算機技術的快速發展，圓柱繞流的研究開始向數值模擬計算方向發展。Frohlich[22]分別采用大渦模擬方法中Smagorinsky和動態Smagorinsky亞網格模型進行數值模擬，Lee等[23]采用有限單元法對相同的模型展開了數值研究，他在圓柱周圍進行了網格加密，得到了與Frohlich相似的流動特征。這些學者們的數值模擬結果雖然得到了與前人試驗相似的流動現象，但是由于計算能力依舊只是計算了無耦合情況下的的流場流動特性。國內的高健停等[24]對有限長三維柱體的截面形狀對其繞流流動特性進行了分析。結果發現自由端的下洗流動對柱體兩側旋渦有較大的抑制作用，這使得阻力系數、升力系數和斯托哈爾數等系數均有較大的減小，而回流區長度則因為受到下洗流動的壓迫而變長，旋渦脫落的模式也有一定的改變，但是這些仿真的研究成果都是重點關注有限長圓柱繞流時的流場特性，沒有考慮結構振動響應和流場耦合的效應。

綜上所述，目前通過試驗的方法，雖然已經成功探索出在高風速下亦可能發生渦激振動，但是并沒有得出確切的高風速下的渦激振動現象發生的機理。在通過數值模擬方法探索的研究領域，由于需要消耗巨大的計算資源和時間，大多數學者僅進行二維的雙向流固耦合分析；針對有限長懸臂圓柱結構，僅計算三維的無耦合固定圓柱的流場流動情況。本文通過商業軟件建立雙向流固耦合數值模擬方法，針對一端固支的圓柱結構在風載荷作用下的頻率響應特性進行分析。通過分析結構振動響應以及漩渦泄落頻率，并與文獻試驗結果[25]對比，表明了用數值模擬來研究普通渦激振動的可行性。通過對圓柱展向各截面上升力曲線進行頻譜分析，發現有限長懸臂圓柱在發生渦激振動時的三維流動特性，揭示了在高風速下也會產生渦激振動的機理。

1 數值模型及求解方法

本文所模擬的風洞試驗模型[25]如圖1所示。在風洞內部受氣流影響的部分為一圓柱，直徑D=0.02 m，高度為25D。圓柱通過風洞地板下的板彈簧固定，能夠產生橫風向的振動，圓柱固有頻率為17.5 Hz，阻尼比為0.28 %。在試驗中監控了在圓柱后方距離圓柱中心5D，偏離圓柱軸線D的沿高度方向上的一系列點的風速波動，從而得到圓柱在不同高度上的渦泄頻率。

圖1 風洞試驗模型[25]Fig.1 Experimental model in wind tunnel experiment

基于計算流體動力學和結構動力學理論，本文在Workbench平臺下，采用ANSYS Mechanical APDL作為結構求解器，CFX作為CFD求解器，圓柱和流體接觸面為耦合面，流場和結構的載荷信息交換發生在流固耦合面上，耦合的邊界發生運動會導致流場網格重構，從而重新進行迭代流場的壓力分布。在求解的每一個分析時間步中，兩個物理場不斷的進行迭代修正，直到在該時間步內滿足收斂條件才會進入下一個時間步。計算開始時，先進行穩態的CFD分析，得到圓柱流固耦合界面的壓力分布，作為圓柱結構瞬態動力學的初始載荷條件，同時穩態流場分布也作為CFD瞬態分析的初始條件。

1.1 流場分析模型

根據風洞試驗模型的尺寸建立流場模型，圓柱中心距離流場進口的長度為10D，距離出口的長度為25D。流場的高度為40D，寬度為20D。整體的流固耦合分析幾何模型如圖2所示。

圖2 流固耦合分析幾何模型Fig.2 FSI analysis geometric model

流體模型的氣動網格是在ANSYS ICEM中劃分，圓柱周圍局部空間的網格采用O-BLOCK方法來進行加密，在保證壁面網格無量綱第一層厚度值y+=u*y/v小于2的前提下，其中，u*為來流速度，y是距離壁面第一層網格高度，v為流體的運動黏度，經過網格無關性測試，最終將壁面第一層網格厚度設為0.006D，以此更加精確的模擬附面層效應。氣動模型網格如圖3所示。

圖3 氣動模型網格Fig.3 Aerodynamic model grids

流體介質為空氣，溫度為25 ℃，密度ρ=1.185 kg/m3，動力黏性系數μ=1.84×10-5kg/ms。進口的邊界條件為Inlet，采取速度進口條件。出口的邊界條件為Outlet，采用壓力出口條件，相對壓力為0 Pa。左右側面以及頂部為Wall，采取的是自由滑移邊界條件。底部的邊界的條件也為Wall，但采用的是無滑移的邊界條件。考慮流場底部的邊界層效應，進口風速在0～200 mm之間隨著高度H的變化而變化，越靠近地面越小，在高度方向200 mm以上采取的是均勻來流速度，如圖4(a)所示。同時，進口湍流度越靠近地面越大，出口湍流度和進口相同。當進口風速為2 m/s時，進口湍流度的分布[25]如圖4(b)所示。

圖4 進口風速及湍流度輪廓線示意圖Fig.4 Profiles of inlet wind speed and turbulence

本文假定流場為不可壓，運動黏度μ和密度ρ均為常數，選取RANS方程作為流體動力學的控制方程，并結合SSTk-ω湍流模型對流場進行求解，采用有限體積法對控制方程進行離散，壓力速度的耦合方式采用SIMPLE算法，對流項為高精度，空間差分格式為二階背風格式，時間差分格式為隱式差分格式。流體動力學控制方程[26]如下

(1)

SSTk-ω湍流模型包括湍動能k和耗散率ω的輸運方程。湍動能k的輸運方程為

耗散率ω的輸運方程為

(4)

φ=F1φ1+(1-F1)φ2

(5)

1.2 結構分析模型

與流場匹配的有限元模型如圖5所示，圓柱底部是采用板彈簧支承，頂端為自由端。彈簧底部固定，另一端與圓柱固接，約束圓柱僅在橫向運動。由于文獻[25]中用于風洞試驗的結構參數沒有給出，只給出了結構的固有頻率fn=17.5 Hz，阻尼比ξ=0.28 %。所以在建模時圓柱的材料剛度接近剛性，調整壁厚、材料密度和板彈簧的剛度，使模型的頻率與風洞試驗模型吻合，并設定同樣的阻尼比。這樣只能保證振動頻率吻合，而振動幅值上有差異。

圖5 結構有限元模型Fig.5 Structural finite element model

本文采用的是4節點的殼單元建模。結構的有限元動力學方程為

(6)

其中，字母上面的點表示對時間的導數，y為立管節點處的位移向量，每個節點處的位移包括該節點處的橫向位移和轉角，{Fy(t)}為立管單元節點處的橫向渦激載荷向量，[M]、[C]和[K]分別為質量矩陣，阻尼矩陣和剛度矩陣，它們分別由各自的單元矩陣集成而得。模型的阻尼采用Rayleigh阻尼，即質量矩陣和剛度矩陣的線性組合

[C]=α[M]+β[K]

(7)

式中：α和β為常數，本文中：α=5.18×10-5，β=4×10-9。

在ANSYS Mechanical APDL中，圓柱橫向運動方程采用顯式動力分析方法進行計算。為保證計算的穩定與結果的準確，還要兼顧計算效率，需要設置合理的時間步長，時間步長是參考庫朗數Co來設定的，同時考慮結構的振動頻率，時間步長設定為0.000 8 s。庫朗數表達式為

(8)

式中：δt為時間步長；U為進口速度；δx為在速度方向的最小網格尺寸。

2 數值計算結果及分析

2.1 結構振動響應分析

圖6是在文獻[25]風洞測得的一系列不同進口風速的工況下，圓柱的無量綱橫向振幅均方根值y/D和約化風速Ur=U/fnD之間的關系。當U=2 m/s時，根據斯特勞哈爾數公式St=fsD/U，可以預估渦泄頻率在17～20 Hz之間，此時約化風速Ur為5.7，圓柱的振幅均方根值出現圖6中第一個峰值，圓柱發生普通渦激振動。當U=5.95 m/s時，預估渦泄頻率在48～57 Hz，Ur為17，此時在曲線上也可見一個局部峰值，圓柱發生高風速下的渦激振動。因此，下面主要針對這兩種進口風速工況下渦激振動響應進行數值模擬研究，通過分析渦泄頻率和結構振動響應頻率之間的關系，探究在高風速下也會發生渦激振動的機理，并與文獻實驗結果[25]進行對比。

圖6 折減振幅(RMS)y/D和進口約化風速Ur的關系Fig.6 y/D versus Ur

通過在Workbench平臺下，針對懸臂圓柱進行雙向流固耦合分析，監控點為圓柱頂部面中心點，當約化風速Ur=5.7時，其位移時間歷程曲線和整個耦合面(即圓柱表面)的升力時間歷程曲線的數值模擬結果如圖7。其中位移為無量綱位移A*=y/D，即真實位移除以直徑D。對圖7的時程曲線進行功率譜密度分析，得到振動響應頻率以及渦泄頻率，如圖8所示。另外，在計算中還發現，當所給風速條件對應的渦泄頻率與結構固有頻率相差稍多時，結構響應會出現明顯的“拍”現象，而圖7(a)中并沒有出現“拍”現象，這表明了該計算成功的捕捉到頻率“鎖定”現象。從圖8(a)可以看出振動響應頻率單一且接近結構的固有頻率17.5 Hz。故認為當Ur=5.7時，結構發生普通渦激振動。

圖7 結構振動響應時程曲線 (Ur=5.7)Fig.7 Time history of vibration response (Ur=5.7)

圖8 功率譜密度分析結果 (Ur=5.7)Fig.8 The result of PSD analysis (Ur=5.7)

另外，從圖8(b)中升力的頻率響應來看，渦泄頻率的大小和圓柱結構的固有頻率相接近，因此，普通渦激振動的發生是由于渦泄頻率靠近圓柱結構的固有頻率，導致了結構發生普通渦激振動。這與文獻中的實驗結果也是十分吻合的。在圖8(b)中，升力的功率譜密度在11 Hz左右有一個小的峰值，這是由于圓柱的橫向振動導致風速的波動而產生的。

當約化風速為Ur=17時，頂部面中心點的位移時間歷程曲線和耦合面上的升力時間歷程曲線的數值模擬結果如圖9所示。對應的功率譜密度如圖10所示。

圖9 結構振動響應時程曲線(Ur=17)Fig.9 Time history of vibration response (Ur=17)

從圖10(a)可以看出，位移振動響應中出現了17.5 Hz和54 Hz，而在耦合面的升力功率譜密度中存在54 Hz和微弱的18 Hz。由圖6中的實驗結果可知，在該風速下也會發生顯著的振動。但是，在圖10 (b)中整個耦合面升力功率譜密度中并沒有出現固有頻率成分。所以進一步分析圓柱高度方向上的渦泄特性。

圖10 功率譜密度分析結果 (Ur=17)Fig.10 The result of PSD analysis (Ur=17)

2.2 流場三維特性分析

為了進一步研究導致高風速下圓柱發生渦激振動的機理，在圓柱展向上取不同的截面進行分析，得到各截面的渦量云圖和升力時間歷程曲線，如圖11～14所示。再針對不同截面的升力時程曲線進行頻譜分析，得到各截面處的渦泄頻率，如圖15所示。

圖11 渦量云圖 (Ur=5.7)Fig.11 Contours of vorticity (Ur=5.7)

圖12 渦量等值面圖(Ur=5.7)Fig.12 Q Iso-surfaces of vortex streets(Ur=5.7)

圖13 渦量云圖 (Ur=17)Fig.13 Contours of vorticity(Ur=17)

圖14 渦量等值面圖(Ur=17)Fig.14 Q Iso-surfaces of vortex streets(Ur=17)

圖15 展向上各截面的渦泄頻率Fig.15 The vortex shedding frequency of each section in the span

從各截面的渦量云圖可以看出，流體流過圓柱時在不同高度處的渦量存在明顯的不同。從基于Q準則畫的渦量等值面圖可以看出，圓柱下游的尾流漩渦脫落形成的渦管沿柱體展向呈現扭曲現象。這都表明了流體流動具有明顯的三維特性。同時，在空氣與圓柱發生流固耦合的過程中，圓柱的振動變形對其尾渦動力特性會產生明顯的影響，不僅會使其尾渦強度發生明顯變化，而且還會使其尾渦分離點發生改變。

圓柱頂端的流線圖如圖16所示，可見氣流呈現明顯的下洗運動，在圓柱頂部產生一個頂渦。由于下洗運動導致了在圓柱頂部區域的漩渦脫落強度明顯減弱，從而在z/L=0.8～1.0處，存在多頻渦泄現象，渦泄頻率中存在18 Hz和54 Hz兩種頻率，18 Hz的頻率和圓柱的固有頻率相近。圖15(d)是z/L=0.8截面處的渦泄頻率。因此可以認為由于流動的三維特性使得接近圓柱自由端處的渦泄頻率接近于結構固有頻率，在高風速下也引起渦激振動，雖然其位移振動幅值小于普通渦激振動的位移幅值，但在工程中，該現象也應該引起重視。

圖16 流場中面上的速度流線圖 (Ur=17)Fig.16 Velocity streamline diagram in the middle of the flow field (Ur=17)

3 結 論

通過數值模擬方法針對有限長懸臂圓柱發生普通渦激振動和高風速下渦激振動的振動機理進行了研究，得到了如下結論：

(1)基于商業軟件ANSYS+CFX建立的數值模擬方法進行雙向流固耦合分析是具有一定的可行性。針對懸臂圓柱進行渦激振動數值模擬結果與文獻實驗結果吻合度較好。

(2)在超過鎖振范圍的高風速下，圓柱結構也會發生渦激振動。雖然其振動幅值小于普通渦激振動引起的振動幅值，但是局部的振幅急劇增加也應該引起重視。

(3)在高風速下的發生渦激振動的原因是由于在結構展向上的三維流動特性引發的，隨著高度的增加，渦泄頻率降低，在z/L=0.8～1.0處，存在多頻渦泄現象，渦泄頻率中低頻成分存在與圓柱結構固有頻率吻合，故此導致了在高風速下也會引起渦激振動。