非對稱粗糙帶對單圓柱流致振動特性影響研究

章大海， 王文顥， 李天嬌， 馮 蕾， 孫 海

（1. 中國石油大學（華東） 化學工程學院， 山東 青島266580； 2. 哈爾濱工程大學航天與建筑工程學院， 哈爾濱150001）

0 引 言

當一定速度的流體繞流柱體結構時， 交替脫落的不對稱旋渦對柱體橫向及流向產生周期性作用力，進而引發柱體振動，柱體振動又反過來影響其尾渦脫落及形態，這種結構與流體之間的相互作用問題就稱為流致振動（Flow-Induced Vibration，FIV）。它是一種廣泛存在于工程應用和自然界中的重要現象，常見的流致振動現象包括渦致振動（Vortex-Induced Vibration，VIV）、顫振（Flutter）、抖振（Buffeting）和馳振（Galloping）等[1]。流致振動會導致高煙囪、海洋油氣輸送管線和海上平臺等建筑物產生疲勞破壞。 在近幾十年核工程、航空航天工程和海洋工程等迅速發展的背景下，流致振動因重大危害性而受到了高度關注[2]，比如1940 年Tacoma 吊橋在較低風速下產生顫振失穩而倒塌。 早在1968 年，Feng[3]就彈性安裝的光滑圓柱渦致振動進行了一些重要的經典測量，得到的振幅比曲線包括“初始分支”和“下分支”，這個定義是由Khalak 和Williamson[4]給出的，他們得到的渦致振動振幅比曲線包括“初始分支”、“上分支”、“下分支”和“去同步化”四部分，得到這種結果的主要原因是Feng 的實驗介質為空氣，而Khalak 和Williamson 的實驗介質為水。 渦致振動可在一定條件下轉化為馳振，馳振有大振幅和小頻率的特點[5]。 Bokaian 和Geoola[6]指出鈍體為非圓柱結構即為方柱等結構可以激發馳振，同時Bokaian和Geoola[7]研究表明當一個圓柱靠近其他圓柱時也可以激發馳振。 Achenbach[8]則表示圓柱外表附有粘貼物可使渦致振動進一步激發為馳振， 這種通過在圓柱表面粘貼附屬物來控制流致振動的方法稱為被動湍流控制（Passive Turbulence Control，PTC）。 為了更好地理解存在PTC 的物理現象，一些研究者把注意力轉向了二維幾何分布，研究帶有對稱或不對稱PTC 的單個圓柱，是了解更復雜的結構所涉及物理現象的一個重要過程。 Nebres 與Batill[9]研究了帶有單個金屬線圓柱的頻率響應，發現旋渦脫落頻率在一定范圍內隨著金屬線分布角度而改變。 Ekmekci 與Rockwell[10]確定了基于近尾特性的單個金屬線圓柱兩個重要分布位置（θc1=55°和θc2=65°），他們發現在Re=10 000 時，金屬線（d/D=0.029）置于θc1可以擴展順流向的尾流氣泡長度，而置于θc2會收縮；他們還發現當金屬線直徑大于未擾動的邊界層厚度時，放置在θc1會明顯減小馮·卡門頻率下的速度波動的頻譜幅度，并且當同一導線處于θc2時，會產生相反的效果。

流致振動不只是會帶來前文提到的危害，它也能造福于人類。 在2006 年，密歇根大學海洋可再生能源實驗室（Marine Renewable Energy Laboratory，MRELab）的Bernitsas 等[11]提出渦致振動水生清潔能源系統（Vortex Induced Vibration for Aquatic Clean Energy, VIVACE），把渦致振動和海洋能發電技術相融合，使這類潛在的危害性轉變為能量的利用。 在之后的研究中，Bernitsas 和Raghavan[12]以選擇性表面粗糙度的方法來加強圓柱的振動，進而增大能量輸出，并且取得了良好的效果。Chang 等[13]對對稱粗糙帶單圓柱的流致振動做了大量研究， 得到了不同規格對稱粗糙帶下的振幅及功率等數據。 2011年，Park[14]等在MRELab 對對稱粗糙帶的分布位置進行大量實驗，進而建立了PTC-TO-FIV 圖，為對稱粗糙帶的選擇提供了參考依據。 近幾十年來，隨著計算機技術飛速發展，計算流體力學（CFD）也被廣泛應用，因為實驗可測量參數有限，CFD 作為補充就必不可少。 Kinaci[15]以k-omega SST 湍流模型研究了光滑單圓柱的流致振動現象，得到了其振幅及頻率曲線，并與Khalak 和Williamson[16]的實驗數據進行了對比，捕捉到了不同的振動區間。 Ding[5]等采用CFD 和實驗研究了對稱粗糙帶單圓柱流致振動特性，對比了兩者的振幅、頻率和功率等，也捕捉到了非對稱的尾渦結構。

然而近年來相關的研究主要集中在雷諾數40 000 以下[17]，更高雷諾數下的流致振動研究較少，對于非對稱PTC 圓柱的研究多見于靜止圓柱，而非對稱PTC 圓柱的流致振動研究較少，在MRELab 進行的圓柱流致振動研究全部集中于對稱粗糙帶，對于非對稱粗糙帶更是尚未見報道。 這種方式是否會造成大的振幅進而產生更大的輸出功率，其尾渦模式又如何，這些都是未知的，而且進一步將對稱和非對稱粗糙帶圓柱及光滑圓柱的流致振動數據進行對比，對更好理解粗糙帶的作用機理有很大幫助，也為理解更復雜結構的流致振動機理做鋪墊。 所以， 本文采用二維非定常RANS 方法和Spalart-Allmaras 湍流模型，結合ANSYS-Fluent 軟件動網格技術和用戶自定義（UDF）功能，對大范圍雷諾數的非對稱粗糙帶單圓柱流致振動進行數值模擬研究， 并與來自美國密歇根大學海洋可再生能源實驗室的實驗結果對比，研究結果具有較為重要的科研意義和廣闊的應用前景。

1 計算模型

本文采用商業CFD 軟件ANSYS-Fluent 對非對稱粗糙帶單圓柱進行數值計算。

1.1 幾何模型

本文的計算區域為矩形，如圖1 所示。 與實驗參數一致[18]，取圓柱直徑D 為88.9 mm，整個計算區域的長為40D=3 556 mm，深為14D=1 244.6 mm，圓柱中心位于入口處15D=1 333.5 mm，并且關于上下壁面對稱。

圓柱的非對稱粗糙帶布置如圖2 所示，圓柱上用高為0.847 mm 的小凸臺來模擬粗糙帶，小凸臺布置位置為α=20°，且覆蓋范圍為θ=16°。

圖1 計算域示圖Fig.1 Computational domain diagram

圖2 粗糙帶分布圖Fig.2 ONE side PTC distribution

1.2 數學模型

對整個計算域二維非定常不可壓縮粘性流體，其雷諾時均方程采用Spalart-Allmaras 湍流模型[19]，可表示為如下方程:

式中:xi為笛卡爾直角坐標，ui為對應的速度分量，ν 為分子運動粘度，t 為時間，ρ 為流體的密度，τij為雷諾應力張量，Sij為平均應變率張量，其表達式為:

1.3 物理模型

本文采用經典的彈簧-質量-阻尼振子模型[20]來描述單圓柱振子的動力特性，單自由度運動方程可以表示為:

式中:m 為系統質量，c 為系統阻尼，k 為系統剛度，FL為圓柱總升力。

一個周期內的平均功率可表示為:

將經典的彈簧-質量-阻尼振子模型的公式（6）代入可得:

為了簡化公式，從理論上求得可獲得的功率，對剛性圓柱的響應假定為近似正弦，則可以表示為:

式中:A 是振幅，δ 是位移曲線的平衡位置與y=0 的偏差，ω 為振動圓柱的角頻率且ω=2πf， f 為圓柱振動頻率。 將位移求導可得速度為:

于是功率可簡化為:

經積分可得VIVACE 系統的輸出功率表達式為:

1.4 邊界條件設定

將所建立的幾何模型劃分為五個區域，如圖3 所示。

圖3 幾何模型整體分塊圖Fig.3 Block diagram of geometric model

圖4 網格示意圖Fig.4 Grid diagram

對如圖3 所示的邊界條件設置為:

（1） 入口為速度入口（Velocity inlet）；

（2） 出口為出流（Outflow）；

（3） 上壁面、下壁面和圓柱為固體面（Wall）；

（4） 圓柱上下位置的實線為內部面（Interior）；

（5） 虛線為交界面（Interface）。

1.5 網格劃分及驗證

本文在整個區域采用二維結構化網格，如圖4 所示，在圓柱四周劃分了2D×2D 的正方形網格區域，將正方形區域使用實面拼接分塊并對圓柱附近進行局部加密，在小凸臺高度方向劃分了兩層。 參考Ding[5]的網格無關性，考慮了相同三種不同疏密程度的網格，對靜止圓柱做了雷諾數為30 000 的網格無關性驗證，其結果如表1。結果表明，三種網格密度的升力系數最大值Cl 和阻力系數平均值Cd 相差較少，為加快計算速度，本文最終采用較稀網格進行計算，網格總數為94 932，Ding 的網格總數為60 252，可認為結果具有網格無關性。

表1 網格無關性驗證Tab.1 Grid resolution study

1.6 UDF 及動網格設定

本文在UDF 中編寫振動方程，對公式（6）采用迭代方法如下:

其中在UDF 里直接賦予m=7.826 kg，k=600 N/m，c=2.64 Ns/m 和c=18.5 Ns/m，FL利用UDF 程序中begin_f_loop( f, t)宏對圓柱表面壓強積分得到，y 利用F_CENTROID(x, f, t)宏函數得到，y˙n和y¨n由公式（13）-（15）得到，在起始時刻定義為y1＝0 及y˙1=0。

導入并加載UDF 之后，開啟動網格，選擇鋪層形式，實現了結構化網格鋪層形式，相對于非結構化其他形式，生成網格簡單迅速，避免了大面積的網格重構，節省了資源的同時又有較高的精度。 其網格運動如圖5 所示，圖5（a）為PTC 圓柱初始時刻的網格圖，圖5（b）為PTC 圓柱運動到某一時刻的網格，可明顯看出運動后網格鋪層產生的網格層。

圖5 網格運動示意圖Fig.5 Grid motion diagram

1.7 參數說明

本文在數據處理中，對一些參數進行了如下無量綱化:

式中:U 為來流速度， fn,w為圓柱固有頻率。

2 結果與討論

2.1 振幅響應

如圖6 為c*=0.02 且U*=10.28 的位移比響應過程曲線。 從圖中可以看出，非對稱粗糙帶圓柱位移比響應過程曲線相對于平衡位置有約0.2 的偏移量，這就表示非對稱粗糙帶圓柱在正負方向的振幅值不同，而且表現為粗糙帶側的振幅小于另外一側，這一結論與對稱粗糙帶圓柱的對稱振幅是不同的。 在本文的數據處理中，為方便比較及討論，后文出現的非對稱粗糙帶振幅均取正負向幅值的平均值。

如圖7 所示為c*=0.02 時不同來流速度的光滑圓柱、非對稱粗糙帶圓柱和對稱粗糙帶圓柱流致振動實驗振幅比曲線，圖中縱軸給出的為振幅比A/D，橫軸為雷諾數Re、來流速度U 和約化速度U*。

圖6 位移比響應過程曲線Fig.6 Displacement response process

隨來流速度增大，光滑圓柱流致振動振幅比曲線依次捕捉到初始分支、上分支、下分支；與光滑圓柱不同， 對稱和非對稱粗糙帶圓柱依次捕捉到初始分支、上分支、過渡分支和馳振分支；粗糙帶對于上分支有明顯的拓寬作用， 光滑圓柱上分支范圍約為5.5≤U*≤7.2，而粗糙帶圓柱約為5.0≤U*≤8.6；與光滑圓柱不同，粗糙帶圓柱振幅比曲線在上分支之后進一步增大，而非進入下分支；在上分支， 非對稱粗糙帶圓柱振幅比大于對稱粗糙帶圓柱， 而在馳振區情況是相反的； 相對于對稱粗糙帶， 可看出非對稱粗糙帶對于分支的分布范圍影響不大。

圖7 振幅響應曲線Fig.7 Amplitude response

進一步對比不同阻尼下的非對稱粗糙帶單圓柱流致振動數值模擬和實驗結果，圖8 為c*=0.02和0.14 時的非對稱粗糙帶單圓柱實驗與數值模擬振幅比曲線。 從圖中可觀察到較明顯的幾部分區域:

（a） 在Re＜40 000～44 000 范圍內捕捉到了初始分支，c*=0.02 時，實驗圓柱的振幅增大到1.1D；c*=0.14 時，實驗圓柱的振幅增大到0.9D。數值模擬比實驗更早地開始振動，主要是由于數值模擬采用簡化的彈簧-質量-阻尼模型，而實驗中的模型為一種更復雜的粘性非線性模型。

（b） 在40 000～44 000≤Re＜70 000～74 000 范圍內捕捉到了上分支，也就是渦致振動（VIV）區，c*=0.02 時，實驗圓柱的振幅保持在1.1D 左右；c*=0.14 時，實驗圓柱的振幅保持在0.9D 左右，而數值模擬的振幅保持緩慢連續增長。

（c） 在70 000～74 000≤Re＜90 000 范圍內捕捉到了上分支-馳振過渡分支，c*=0.02 時，隨著Re的增大，實驗圓柱的振幅比從1.1D 增大到1.2D；c*=0.14時，實驗圓柱的振幅0.8D 減小到0.7D。 圖中給出的為多次實驗的平均結果， 在此區間內振幅比標準偏差相對較大，說明過渡期振動不穩定，這就造成數值模擬捕捉的困難， 數值模擬在此區間和實驗吻合不好也能歸因于URANS 湍流模型不能預測隨機小渦。 需要指出的是，在Kinaci[21]和Wu[22]等模擬研究中，過渡區也沒有很好地吻合。

（d） 在90 000≤Re＜120 000 范圍內捕捉到了馳振分支， 隨著Re 的增大，c*=0.02 時圓柱振幅從1.2D迅速增大到1.7D；c*=0.14 時從0.7D 迅速增大到1.3D。在馳振區前段，實驗與數值模擬吻合良好，后段實驗圓柱的振幅小于數值模擬的振幅是因為實驗中到達了安全停止位置，這樣做的目的是保護實驗裝置，同時實驗中自由液面的影響也限制了圓柱振幅的進一步增加。

圖8 振幅比曲線Fig.8 Amplitude ratio curve

從整體來看，數值模擬結果和實驗結果大體吻合，能夠很好地表現出與實驗相對應的整個來流速度變化范圍內的振幅變化，在每種系統阻尼比下，非對稱粗糙帶單圓柱流致振動實驗振幅比曲線依次經歷初始分支、上分支、上分支-馳振過渡分支和馳振分支，但數值模擬捕捉上分支不明顯；阻尼增大使得非對稱粗糙帶單圓柱流致振動的振幅下降。 另外，在進入馳振分支后的振幅突然增大也符合馳振大振幅的特點。

2.2 頻率響應

圖9 頻率比曲線Fig.9 Frequency ratio curve

圖9 為c*=0.02 和0.14 時的實驗與數值模擬頻率比曲線。 頻率比曲線的區域與振幅比曲線對應，數值模擬結果和實驗結果趨勢大體吻合，除上分支區外，頻率比捕捉較為準確；阻尼對于非對稱粗糙帶單圓柱流致振動頻率的影響不大。在上分支區，實驗與數值模擬存在較大誤差，在眾多的研究中，這是雷諾時均方法存在的缺陷，這也是在相關領域模擬中廣泛存在的問題，在Kinaci[21]和Wu[22]等文章中存在同樣的問題，圖9 中給出了Kinaci 的數值模擬曲線。 另外，在進入馳振區，實驗中的f*下降并穩定于1.1，數值模擬的f*下降并穩定于1，符合馳振小頻率的特點。

2.3 尾渦模式

關于尾渦模式的定義最早由Williamson 和Roshko[23]給出，包括2S、2P 和2C 等，S 表示單個旋渦，P 表示一對旋轉方向相反的旋渦，C 表示一對旋轉方向相同的旋渦。 研究表明尾渦模式主要取決于雷諾數，如圖10（a）-（f）為c*=0.02 時不同雷諾數下的尾渦流態圖，其中紅色和藍色表示旋渦方向相反。

隨著雷諾數的增大，圖10(a)所示的2S 模式逐步過渡為如圖10（f）所示的2P 模式，尾渦被逐漸拉長，旋渦脫落距離圓柱尾部越遠，拉長的旋渦具有更多的能量，因而可以激發圓柱振子更大的振幅。

取c*=0.02，U*=13.37，對應的Re =124 460，這時已處于馳振區，對一個周期內的各個時刻尾渦流態圖進行研究。圖11（a）-（g）為同一周期內不同時刻的尾渦流態圖。可以發現的是，一個周期內圓柱尾部共有兩對相反的旋渦脫落，屬于2P 型，但與傳統2P 型不同的是，在上側脫落的旋渦對大小基本相等，如圖11（d）黑色線框所示，而下側脫落的旋渦對大小不等，如11（b）黑色線框所示。

圖10 c*=0.02 的尾渦形態Fig.10 Wake vortex (c*=0.02)

圖11 c*=0.02 的馳振區尾渦形態Fig.11 Wake vortex of galloping (c*=0.02)

2.4 功率計算

利用前文給出的功率計算公式得出功率， 進而畫出曲線， 圖12 所示為c*=0.14 的功率曲線。 c*=0.02 時，所得功率用于系統阻尼耗散，并不能得出輸出功率，所以未畫出。

由圖可看出，同對稱粗糙帶圓柱相比，非對稱粗糙帶圓柱流致振動在上分支有更大的輸出功率，在馳振分支則結論相反。對稱粗糙帶c*=0.14 時，非對稱粗糙帶單圓柱最大輸出功率為6.62 W，實驗所得最大功率為7.29 W。 除上部分支外，趨勢吻合良好，由前文所給計算公式及第2.2 節頻率響應可看出上分支區吻合不好主要可歸因于頻率捕捉的不準確。

3 結 論

本文利用ANSYS-Fluent 軟件模擬非對稱粗糙帶單圓柱流致振動現象， 結合UDF 功能和動網格技術來控制非對稱粗糙帶圓柱的振動，研究討論了不同阻尼下的振幅、頻率、尾渦模式和功率輸出等，并與實驗結果進行對比。 得出的結論主要有:

（1） 非對稱粗糙帶單圓柱流致振動正負向振幅值不同，粗糙帶側位移值較小。

（2） 非對稱粗糙帶振幅比曲線捕捉到了初始分支區、上分支、上分支-馳振過渡分支和馳振分支，但數值計算捕捉上分支不明顯。

圖12 c*=0.14 的功率曲線Fig.12 Power curve (c*=0.14)

（3） 相對于光滑圓柱，非對稱粗糙帶對上分支有明顯的拓寬作用。

（4） 與光滑圓柱不同，非對稱粗糙帶單圓柱流致振動在進入馳振分支后，頻率比下降并穩定于某一值，本文實驗為1.1，數值計算為1。

（5） 在所研究范圍內，隨著雷諾數的增大，非對稱粗糙帶單圓柱的尾渦被逐漸拉長，旋渦脫落距離圓柱尾部越遠，模式從2S 過渡為非對稱的2P。

（6） 同對稱粗糙帶圓柱相比，非對稱粗糙帶圓柱流致振動在上分支有更大的輸出功率，馳振區則相反。

（7） 本文通過數值計算得到的非對稱粗糙帶單圓柱最大輸出功率為6.62 W，實驗所得最大輸出功率為7.29 W。