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

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

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

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

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

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

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

1 數(shù)值模型及求解方法

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

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

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

1.1 流場分析模型

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

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

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

圖3 氣動模型網(wǎng)格Fig.3 Aerodynamic model grids

流體介質(zhì)為空氣，溫度為25 ℃，密度ρ=1.185 kg/m3，動力黏性系數(shù)μ=1.84×10-5kg/ms。進口的邊界條件為Inlet，采取速度進口條件。出口的邊界條件為Outlet，采用壓力出口條件，相對壓力為0 Pa。左右側(cè)面以及頂部為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

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

(1)

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

耗散率ω的輸運方程為

(4)

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

(5)

1.2 結(jié)構(gòu)分析模型

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

圖5 結(jié)構(gòu)有限元模型Fig.5 Structural finite element model

本文采用的是4節(jié)點的殼單元建模。結(jié)構(gòu)的有限元動力學方程為

(6)

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

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

(7)

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

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

(8)

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

2 數(shù)值計算結(jié)果及分析

2.1 結(jié)構(gòu)振動響應分析

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

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

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

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

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

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

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

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

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

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

2.2 流場三維特性分析

為了進一步研究導致高風速下圓柱發(fā)生渦激振動的機理，在圓柱展向上取不同的截面進行分析，得到各截面的渦量云圖和升力時間歷程曲線，如圖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準則畫的渦量等值面圖可以看出，圓柱下游的尾流漩渦脫落形成的渦管沿柱體展向呈現(xiàn)扭曲現(xiàn)象。這都表明了流體流動具有明顯的三維特性。同時，在空氣與圓柱發(fā)生流固耦合的過程中，圓柱的振動變形對其尾渦動力特性會產(chǎn)生明顯的影響，不僅會使其尾渦強度發(fā)生明顯變化，而且還會使其尾渦分離點發(fā)生改變。

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

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

3 結(jié) 論

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

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

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

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