質量比對D 形截面柱體流致振動的影響1)

宋吉寧 李 壯 蔣學煉 金瑞佳 劉宇航

* (天津城建大學天津市軟土特性與工程環境重點實驗室,天津 300384)

? (交通運輸部天津水運工程科學研究院,港口水工建筑技術國家工程研究中心,天津 300456)

引言

流致振動是一種復雜的流固耦合現象[1],廣泛存在于海洋立管、海底管道、大跨橋梁結構等工程領域中,可導致結構疲勞破壞.已有的研究大多關注于認識和抑制結構的流致振動,減少流致振動所帶來的危害[2-4].近年來,隨著綠色能源開發研究的發展,利用流致振動實現海流能俘獲也引起了學者們的關注[5-9],具有代表性的是Bertinsas 等[10-11]開發了一種低速海流發電裝置(VIVACE)進行海流能的捕獲,此裝置可將柱體振動的機械能轉化為電能,具有啟動流速低、簡單易維護、對海洋環境好、不影響通行等潛在優勢.關于海流流致振動和海流能能量轉化的研究,練繼建等[12]梳理了柱體流致振動的研究現狀并剖析其中存在的問題,建議可利用流致振動進行新型能源的開發.

在柱體流致振動和能量捕獲的研究中,大多數是關于圓形截面柱體的[13],及春寧等[14]和殷布澤等[15]對圓柱形的研究進行了綜述.Williamson 等[16-18]開展模型實驗研究了低質量比圓柱的渦激振動特性,發現振幅與整個系統的質量阻尼比(m?ξ)有關.宋吉寧等[19]使用拖曳裝置在均勻流條件下對柔性立管進行渦激振動實驗,發現在均勻流條件下,立管仍存在多模態振動,順流向主導模態則多達12 個.陳正壽等[20]針對質量比開展了剛性圓柱的數值模擬,發現質量比是影響圓柱體振幅的一個重要因素.唐國強等[21]發現在低雷諾數下,當串聯的柱體距壁面的距離/柱徑為0.25 時,柱體的尾渦脫落被抑制,柱體不發生振動.劉旭菲等[22]對不同質量比的圓柱進行了近壁面渦激振動的數值模擬.發現隨著質量比的增大,近壁面圓柱的振動在高折合速度下才會發生,振幅也較小,且受到壁面邊界層的影響,低質量比柱體在順流向和橫流向的振動頻率相等.在能量俘獲效率方面,白旭等[23]發現圓柱體的獲能效率與柱體振幅的大小并沒有直接的關系,柱體的振幅越大,獲能效率不一定越高.

已有研究表明,不同的截面形狀對柱體流致振動有重要影響[24].李海濤等[25]通過實驗和CFD 的方法研究了不同截面下鈍頭體及寬厚比對流致振動能量收集特性的影響,發現鈍頭體為D 形截面時,流致振動呈現“渦激振動”“渦激振動-馳振”和“馳振”的變化.

在非圓形截面中,由一個平面和一個半圓的曲面組成的D 形截面柱體(圖1)是比較有代表性的一種.關于D 形截面柱體流致振動的研究,部分是通過風洞試驗[26],然而由于空氣和水在密度和黏性方面差異較大,柱體在水中的流致振動與風洞實驗中存在較大差異.Zhao 等[27]通過水槽實驗研究了質量比為6 時D 形截面柱體的流致振動響應,D 形的直邊分為迎著或背向來流方向兩種狀態,并指出沒有后體也可以發生渦激振動.Chen 等[28]開展了低雷諾數100 時,質量比為2 時單自由度D 形截面柱體橫流向流致振動的數值模擬,分析了不同迎流攻角、渦激振動與馳振的特性.衛昱含[29]結合高斯過程回歸方法分別對質量比為2 的迎流攻角0°和180°的D 形截面雙柱體進行數值模擬發現,柱體的振幅和能量利用效率無必然聯系,角度和阻尼比對能量利用效率影響顯著.Li 等[30]設計了ODO,ODODO 和DOD 三種不同排列的柱體,研究截面變化對渦激振動能量收集的影響.結果表明,D 形截面柱體夾在兩個圓形截面柱體之間時,柱體的振動受到抑制;而其余兩種排列的柱體可以增強渦激振動,并且增大鎖定區間,提高了能量轉化效率.

圖1 D 形截面柱體布置圖Fig.1 Schematic of D-section prism

上述研究表明,關于不同質量比下D 形截面柱體流致振動響應及其尾流特性的認識尚未完善,諸如響應特征、尾渦脫落形態、能量轉化效率等特性或規律尚需進一步研究.因此,本文通過數值模擬研究四個不同質量比(2,5,7 和10) D 形截面柱體的流致振動問題.在雷諾數范圍為288～2880，約化速度范圍為2～20 的條件下,分析柱體橫流向振動響應、頻率響應、脫渦模式以及能量轉化效率等規律特性,以期為海流能開發利用相關工程應用提供參考數據.

1 數值模型與計算方法

1.1 數值模型

在流場中放置一個D 形截面柱體,如圖1 所示圓弧面向下,直邊與水流方向平行,其直邊長度DL為0.0381 m,迎流角度 α (見圖1)固定為 90?.為了探究柱體在橫流向的最大振幅,忽略阻尼的影響,設阻尼比 ξ=0.選取了四個質量比m?=m/md,分別是2,5,7 和10,研究D 形截面柱體在均勻來流U下的橫流向振動特性.柱體在靜水中的固有頻率設置為fn=0.4 Hz,柱體在迎流面的特征長度De(De=0.5(1+|cosα|)DL) 取迎流面投影最大寬度,為0.01905 m.

計算域的參數設置,如圖2 所示,坐標原點位于D 形截面柱體直邊長度的中點處,x軸正方向為順流方向,y軸則垂直于來流方向.流場計算域的順流向總長度為48DL,橫流向總寬度為24DL,柱體直邊中心距離上下兩側邊界及入口邊界均設置為12DL,阻塞率為2.08%,小于5%可以忽略流場寬度對柱體振動的影響[31].

圖2 計算區域與邊界條件設置Fig.2 Computational domain and boundary conditions

本文所模擬的雷諾數范圍為288 ≤Re≤2880,入口流速為恒定流.出口采用壓力出口的邊界條件,上下兩側邊界設置為對稱邊界條件,柱體表面設置為無滑移壁面.

本文基于二維非定常不可壓縮流體RANS 方程,采用ANSYS FLUENT 軟件,應用k-ω SST 湍流模型,SIMPLEC 方法求解壓力速度耦合方程,壓力項采用PRESTO 方法離散[32],時間項采用QUICK方法離散,湍流動能項求解采用二階迎風格式,比耗散率項求解采用一階迎風格式,瞬態項采用二階隱式求解方法,各個參數的收斂殘差設置為1.0×10-5,最大迭代步數設置為30 步.

1.2 柱體結構運動計算方法

本文僅研究柱體在橫流向上的振動響應,因此,流致振動的系統可以看作是一個單自由度的彈簧-質量-阻尼模型(圖2),其在橫流向的結構動力學方程為

式中,y為柱體中心在橫流向的位置,m,c,k分別為系統的質量、阻尼和彈簧剛度.其中彈簧剛度k=,ω0為系統的圓頻率.阻尼,ξ 為系統的阻尼比.Fy(t) 為柱體在橫流向的阻力和升力.

柱體運動方程求解采用Newmark-β法,計算出每一個時間步后柱體新的位置、速度及加速度.需要注意的是,應用FLUENT 軟件對柱體邊界上的網格節點進行更新時,UDF (用戶自定義函數)設置中不能將柱體新時刻的速度直接傳遞給柱體邊界網格,即vel[1]≠y(t+?t) .否則,將會造成新時刻柱體邊界網格的位置與柱體位置出現偏差[33].根據Newmark-β法已知t時刻柱體的位置,則t+?t時刻柱體的位置為

式中,yt+?t為柱體t+?t時刻位置,yt為柱體t時刻位置.?t為時間步長.根據Newmark-β法無條件穩定的假設,β 一般取0.25,則式(2)中t+?t時刻的柱體的位置為

通過DEFINE_CG_MOTION 宏,可將UDF 中計算出柱體在t+?t時刻的速度,傳回FLUENT 對柱體位置進行更新.根據軟件手冊,FLUENT 使用下式對網格節點位置進行更新

比較式(5)和式(3),發現柱體邊界上網格位置與柱體位置的二階項存在差異,這會導致兩者位移不同步.因此,需要在UDF 中對新時刻柱體邊界網格節點速度計算進行修正,采用式(6)

將上式代入(4),得到柱體邊界網格節點在新時刻t+?t時的位置為

比較式(7)與式(3),每次時間步更新后,柱體邊界網格節點的位移與柱體是一致的.這表明采用式(7)更新柱體邊界網格節點的新時刻速度,可以保證柱體邊界網格與柱體位移同步,以獲得更精確的數據.

1.3 能量轉化效率

在水流作用下D 形截面柱體吸收流體能量產生運動,在一個周期T中,柱體俘能功率可以表達式為[34]

式中,m為柱體的質量,為柱體在橫流向上的速度.流體掃過D 形截面柱體的能量可以表示為

其中,L為柱體的長度,U為來流速度,DL為D 形截面柱體直邊長度.則一級能量轉化效率 η 可以表示為

1.4 網格劃分及網格獨立性驗證

D 形截面柱體的流致振動特性復雜,為了精確捕捉D 形截面柱體周圍及尾部的流動特征,采用結構化網格剖分計算域(圖3),將流體域劃分為9 個部分,D 形截面柱體附近6 倍直徑DL的范圍采用“O 形切分”,近壁面進行局部加密,使得y+<0.6,y是從第一層網格中心到柱體壁面的距離,τw是圓柱壁面的剪切應力,μ 是動力黏性系數,ρ 是流體的密度).

為驗證網格無關性,對比了5 組不同疏密的網格設置.在約化速度Ur=5,m?=2.6,fn=0.4 Hz,ξ=0.003611 時計算D 形截面柱體的橫流向響應振幅和頻率.Mesh1～Mesh5 的主要區別是,逐次減小柱體周圍6DL范圍內的圓形加密區的網格增長率(數值見表1),可使得圓形加密區內的網格數量逐次增大,而離柱體較遠區域的網格數量僅隨之略有增多.總網格數量從2.6 萬增至6.1 萬.

表1 網格獨立性驗證Table 1 Mesh independency study

表1 給出了5 種網格的計算結果,其中Ay/De為柱體無量綱橫流向的均方根振幅(Ay=yrms),fs/fn為柱體振動頻率與固有頻率的比值.通過對比表1 的數據,在網格數量足夠多的情況下,Mesh3,Mesh4 與網格Mesh5 相比,柱體的均方根振幅與振動頻率的變化幅度均在1%以下.兼顧計算精度和效率,本文選擇Mesh3 作為后續的計算網格.

為驗證時間步長 ?t的無關性,在U?t/DL≤0.01的范圍內,對比了5 個不同的時間步長的計算結果,具體參數見表2.算例參數同樣采用約化速度Ur=5,m?=2.6,fn=0.4 Hz,ξ=0.003611 時D 形截面柱體的橫流向響應振幅和頻率結果.從表2 可見,與最小的時間步長工況C1 相比,C4 的振幅偏差0.75%,頻率偏差0.43%,相對偏差幅度都在1%以下,能夠滿足計算精確要求,兼顧計算效率考慮,本文選擇將時間步長設為0.005 s.

2 數值模型驗證

為了驗證本文數值模型的可靠性,采用與D 形截面柱體流致振動實驗相同的計算參數[27],柱體迎流角度為0°,且僅在橫流向上振動,D 形截面柱體的直邊長度DL=0.025 m,柱體的質量m=0.9018 kg,ξ=0.00151,在空氣中的固有頻率為fn,air=0.783 Hz,在水中的固有頻率fn,water=0.74 Hz,附加質量為ma=((fn,air/fn,water)2-1)m.圖4 比較了無量綱振幅(y10/DL)的數值結果和實驗數據,其中y10表示柱體橫流向的前10%最大振幅的平均值.從圖中可以看到數值模擬柱體振動的幅值與實驗值整體吻合較好,均隨著約化速度的增加而增大,驗證了本文數值模型的準確性.

圖4 響應振幅的對比Fig.4 Comparison of the vibration amplitude

3 數值計算結果與分析

3.1 振幅響應

在水流作用下,D 形截面柱體橫流向振動的平衡位置逐漸脫離了柱體靜止時的位置,因此采用柱體振動的平衡位置偏移量[35]和振幅[36]兩者結合進行對比分析.|y0| 為柱體振動時振動平衡位置的偏移量.從圖5 可以看到,當柱體的質量比為m?=2,5時,在低約化速度下(2 .0 ≤Ur≤6.5),柱體處于渦激振動的初始分支和上端分支,此時柱體振動的平衡位置并未出現較大的偏移;隨著約化速度的增大,柱體進入到渦激振動-馳振共同作用的區間,振動的平衡位置逐漸偏離了原來靜止時的位置,向一側偏移.其中質量比m?=2 的柱體,柱體振動的偏移量最大.當柱體的質量為m?=7,10 時,平衡位置偏移量(|y0|/D)隨著約化速度的增加而不斷偏離靜止時的位置,雖然質量比較大,但是柱體偏移量的增長幅度較小.這是由于D 形截面柱體在橫流向上不對稱,升力總是更傾向于指向截面的直邊一側,升力的時間平均值不在零線附近.

圖5 振動平衡位置偏移量隨約化速度的變化Fig.5 Variation of equilibrium position offset with reduced velocity

圖6 為不同質量比下的橫流向D 形截面柱體的最大振幅ymax隨約化速度的變化情況.從圖6 中可知,質量比對柱體振動特性具有顯著的影響.以質量比m?=10 為例,當質量比m?=10 時,可以看到D 形截面柱體振幅響應的四個分支:初始分支(2 ≤Ur≤3)、上端分支(3 ≤Ur≤4.5)、渦激振動-馳振分支(5 .0 ≤Ur≤13.5)以及馳振分支(13.5 ≤Ur≤20).其中柱體的初始、上端分支所在的區間都很窄,柱體在上端分支的最大無量綱振幅只有0.2De,在離開上端分支后,柱體進入了渦激振動-馳振分支,D 形截面柱體與圓柱體的振動特性[28]不同,振幅沒有出現振幅驟降的現象,而是表現為逐漸下降.在Ur=13.5之后,結合柱體振動的頻率圖(圖7),柱體振動的頻率很小(fs/fn,water=0.67),但振幅不斷增大,呈現出高振幅、低頻率特征,此時柱體進入到馳振分支.當柱體的質量比m?=7 時,柱體在Ur=15.0 之后才進入到馳振分支中;當m?=5 時,柱體在更大約化速度時才進入到馳振分支.在本次模擬中,當m?=2 時沒有捕捉到馳振分支,且振幅最小;質量比m?=5 時,柱體振幅出現最大值.從對比來看,在相同的響應分支區間里,質量比越小,振幅往往越大.

圖6 最大振幅隨約化速度的變化Fig.6 Variation of maximum amplitude with reduced velocity

圖7 主導振動頻率隨約化速度的變化Fig.7 Variation of dominant frequency with reduced velocity

3.2 振動頻率

圖7 展示了不同質量比下D 形截面柱體的無量綱振動頻率(fs/fn,water)隨著約化速度的變化趨勢,其中fs為柱體的振動頻率,可通過柱體位移時程曲線進行快速傅里葉變換(FFT)并提取主導頻率后得到.圖8 為D 形截面柱體的在幾個不同約化速度下振動的頻率譜.

圖8 D 形截面柱體振動的頻率譜Fig.8 Vibration frequency spectra of D-section prism

從圖8 可以看到,在低約化速度下,隨著質量比的增大,柱體的無量綱振動頻率隨著約化速度的增加而不斷增大.當柱體的質量比m?=2 時,在約化速度 2 .0 ≤Ur≤4.0 時,柱體處于渦激振動的初始分支和上端分支,柱體的振動頻率是單倍頻且在St=0.2附近,此時柱體的振幅和平衡位置偏移量都很小.在約化速度Ur≥6.5 時,柱體的振動頻率是雙倍頻且在St=0.2 的下方,并且靠近St,振幅很小,近似于渦激振動.但并不能說明此時柱體發生了渦激振動,此時柱體的振幅是不斷變化的,實際上此時的柱體處于渦激振動-馳振共同作用的分支.

當柱體的質量比m?=10 時,柱體在低約化速度下,進入渦激振動分支.在 2 .0 ≤Ur≤3.0 區間中,柱體處于渦激振動的初始分支,柱體的振動頻率落在了St上,此時柱體的振幅很小.在 3 .5 ≤Ur≤6.0 時,柱體進入到了上端分支,振動的頻率比在1.0 附近,結合前文的振幅響應圖可以發現,此時柱體有最大的振幅.此后柱體進入到渦激振動-馳振共同作用的分支,在這個分支中,渦激振動占主導作用,表現為圖8 (m?=10,Ur=11.5)中只有一個主頻.在高約化速度下(Ur≥13.5)進入到馳振分支,此時柱體的無量綱振動頻率均小于1,且柱體的振幅不斷增大,表明柱體進入到了馳振分支.

在同一個質量比的情況下,隨著約化速度的增大,D 形截面柱體經歷了渦激振動初始分支、上端分支、渦激振動-馳振分支以及馳振分支;柱體的振動頻率逐漸由單倍頻向多倍頻轉變;馳振在流致振動中的影響逐漸增大.對比來看,質量比對D 形截面柱體的振動響應有較大影響,質量越大的柱體其渦激振動的初始分支與上端分支的區間越大,柱體越容易進入到馳振分支,但振幅及平衡位置偏移量越小.

3.3 尾渦脫落模式

D 形截面柱體與圓柱相比,柱體的物理幾何尖角對柱體后方尾跡區域的旋渦形成、脫落密切相關.柱體的振動振幅及頻率是由柱體振動時所受到的流體力決定的,同時柱體表面的流動分離剪切層直接影響流體力,柱體后方旋渦的形成與脫落是這兩者共同作用的結果.在本文研究的雷諾數范圍內,根據柱體后方尾渦的形成與脫落過程,觀察到多種不同的旋渦脫落形態,如“2S”(圖9)、“4S+4S”(圖10)等,其中“4S+4S”表示在一個振動周期內柱體兩側各脫落4 個旋渦[28].

圖9 m*=2 柱體的尾渦脫落模式(Ur=13.0)Fig.9 Vortex shedding mode of the cylinder at m*=2 and Ur=13.0

圖10 m*=10 柱體的尾渦脫落模式(Ur=16.0)Fig.10 Vortex shedding mode of a cylinder at m*=10 and U r=16.0

圖9 和圖10 分別給出了m*=2,10 的多個瞬時尾渦脫落的形態.當柱體質量比m?=2,在約化速度Ur=2.0 時(雷諾數為288),D 形截面柱體后方沒有發生明顯的尾渦交替脫落,而在其他的約化速度下,尾渦脫落模式為“2S”模式,如圖9 所示.流體流經柱體時,柱體表面形成上下兩側的分離剪切層,分離剪切層之間相互作用,直邊一側的邊界層上形成的旋渦S1,由于柱體截面存在幾何尖角,在流動的過程中被切斷并脫落,旋渦S1 與S2 大小幾乎相等、旋轉方向相反,交替從柱體上脫落.當柱體質量比m?=10,在約化速度Ur=2.0 時(雷諾數為288),柱體振幅很小,再結合圖7 的振動頻率可見,此時柱體類似于固定圓柱繞流;在約化速度 2 .5 ≤Ur≤5 時,柱體處于渦激振動的分支,尾渦脫落模式為“2S”;在約化速度 5 ≤Ur≤13 時,柱體處于渦激振動-馳振共同作用的區間,且渦激振動占主導;在約化速度13.5 ≤Ur≤20 時,柱體進入到馳振區間,無量綱振動頻率降低至0.67 附近,尾渦脫落模式由“2S”模式轉變成“4S+4S”,即在一個振動周期內從D 形截面柱體兩側各脫落4 個旋渦,如圖10 所示.采用圖9 和10 所示的方法,可以識別各個工況下的尾渦脫落模式,匯總結果見圖11,其中紅色叉號代表著此時柱體無尾渦脫落,“SS”表示柱體兩側后方同時脫落一個旋渦[32].從圖11 可以看到,不同質量比柱體的尾渦脫落模式,在中低約化速度是比較接近的,多數情況是“2 S”模式,但是在高約化速度時,較高質量比的尾渦脫落模式出現了更多的變化,而低質量比m*=2時在本文計算范圍內尾渦脫落形態都表現為“2S”模式.

圖11 4 種質量比不同約化速度的尾渦脫落模式變化Fig.11 Variations of vortex shedding patterns with reduced velocity at four different mass ratios

3.4 一級能量轉化效率

根據前文1.3 部分給出的方法和公式,計算了四個不同質量比下D 形截面柱體流致振動能量的一級轉化效率,結果見圖12.

圖12 能量轉化效率的對比Fig.12 Comparison of the energy conversion efficiency

從圖12 可以發現,在約化速度較低時,由于柱體的振幅、頻率都較低,因此其俘能效率處于較低的水平;而當來流速度增大后,柱體的能量轉化效率也隨之同步增加,柱體的振動頻率逐漸接近系統的固有頻率,所處的區間為渦激振動的區間.隨著約化速度持續增大,馳振對柱體的影響逐漸增強,柱體此時處于渦激振動-馳振的區間,柱體的振幅逐漸降低,能量轉化效率也逐漸下降,當柱體進入到馳振區間后,柱體的能量轉化效率也逐漸增大.出現這種變化的原因,可能是由渦激振動向馳振轉化的過程中,柱體的振動頻率降低導致轉化效率下降,在進入到馳振區間后,D 形截面柱體的振幅發生了大幅增長,從而提高了能量轉化效率,但由于柱體頻率較低,柱體能量轉化效率的提升并不顯著.

通過圖12 的對比來看,質量比對D 形截面柱體的能量轉化效率的影響較為明顯.另外,結合前文3.1 部分關于振動分支的劃分,當柱體處于渦激振動的分支時,其能量轉化效率最大.在本文研究范圍內,質量比m?=10 且Ur=4.5 時,D 形截面柱體一級能量轉化效率達到最大值44%.

4 結論

本文對4 種不同質量比D 形截面柱體的流致振動進行了數值模擬,處于亞雷諾數區間(288 ≤Re≤2880),分別對不同約化速度下的D 形截面柱體的振幅、頻率、平衡位置偏移量、尾渦脫落形態及能量轉化效率等進行了對比分析.本文得到以下結論.

(1)在模擬計算的范圍內,質量比越大的柱體,進入馳振區間對應的約化速度越低.隨著約化速度的變化,D 形截面柱體在低約化速度下均觀察到典型的渦激振動分支,包括渦激振動的初始分支、上端分支.質量比為2 時,隨約化速度增大離開上端分支進入渦激振動-馳振分支,未發現明顯的馳振分支;質量比為5,7 和10 時,柱體離開渦激振動-馳振分支后,進入馳振區間時的約化速度依次減小.

(2)柱體的質量比對柱體平衡位置偏移量有較大的影響.在本文考慮的范圍中,質量比越大的柱體,流致振動的平衡位置偏移量越小.

(3)在渦激振動區間,柱體的尾渦脫落模式多數為“2S”,在馳振區間,尾渦脫落轉為“nS+mS”模式.在本文模擬的工況中,柱體質量越大,尾渦處于“2S”模式的約化速度區間越窄.

(4)約化速度顯著影響柱體的一級能量轉換效率.高能量轉化效率出現在渦激振動分支,而不是在馳振分支.質量比對一級能量轉化效率影響在約化速度2～8 區域較明顯,而在其他約化速度區域的差異較小.

限于篇幅,本文只關注于一個迎流角.至于不同迎流角、不同質量阻尼比的影響以及其他雷諾數范圍的情況,還有待進一步的研究.

致謝

感謝國家超級計算天津中心、天津城建大學高性能計算平臺對本文數值計算工作給予的支持.