康 莊， 張立健， 倪問池

(哈爾濱工程大學(xué) 深海工程技術(shù)中心， 哈爾濱 150001)

在深海工程領(lǐng)域，圓柱體結(jié)構(gòu)例如海洋立管是工程上應(yīng)用最多最廣泛的一種結(jié)構(gòu)物，當(dāng)水流繞過圓柱體時(shí)，會發(fā)生結(jié)構(gòu)物后漩渦脫落的現(xiàn)象，交替的漩渦會激發(fā)與來流方向垂直的流體作用力，致使結(jié)構(gòu)物產(chǎn)生振動，當(dāng)泄渦頻率與固有頻率接近時(shí)，將會發(fā)生“共振”現(xiàn)象，同時(shí)圓柱振動幅值將增大[1]。圓柱體渦激振動(Vortex-Induced Vibration, VIV)一直都是研究所關(guān)注的熱點(diǎn)問題。國內(nèi)外眾多學(xué)者進(jìn)行了彈性支撐圓柱渦激振動研究。黃維平等[2-4]基于流固耦合的渦激升力模型和Morison方程，提出了一個(gè)非線性的圓柱體渦激振動時(shí)域分析模型，并且采用ANSYS-CFX 軟件對細(xì)長柔性圓柱體雙自由度渦激振動進(jìn)行了深入研究。Jauvtis等[5]通過實(shí)驗(yàn)發(fā)現(xiàn)，當(dāng)質(zhì)量比為7時(shí)，單自由度與雙自由度渦激振動試驗(yàn)結(jié)果類似，其橫向最大振幅為1.05D附近，當(dāng)質(zhì)量比為2.6時(shí)，雙自由度橫向振幅出現(xiàn)“超上端分支”，其幅值達(dá)到1.5D，彈性支撐雙自由度渦激振動出現(xiàn)如此大橫向振幅，此現(xiàn)象是以前從未發(fā)現(xiàn)的。谷家揚(yáng)等[6]希望采用數(shù)值模擬方法研究“超上端分支”現(xiàn)象，然而即使將無因次阻尼比ξ設(shè)為0，模擬出的最大橫向振幅也不超過1D。部分學(xué)者注意到, 試驗(yàn)結(jié)果或者計(jì)算結(jié)果跟初始條件有很大的關(guān)系，Guilmineau等[7]參照J(rèn)auvtis等的單自由度低質(zhì)量比試驗(yàn)，計(jì)算了2D橫向彈性支撐圓柱(m*=2.4,m*ξ=0.013)的渦激振動,數(shù)值模擬采用了三種初始條件:①速度保持均勻不變;②速度由最高逐漸減小;③速度由最低逐漸增加。計(jì)算結(jié)果表明：圓柱橫向振幅與初始條件有關(guān):僅在①和②條件下預(yù)測到了“下支”響應(yīng);在條件③時(shí),觀察到了“上支”。試驗(yàn)采用了勻加速和勻減速的形式，從而產(chǎn)生了“遲滯現(xiàn)象”，而數(shù)值模擬一般采用固定速度形式，遲滯現(xiàn)象對渦激振動響應(yīng)有顯著影響。

CFD(Computational Fluid Dynamics)模擬采用固定速度初始條件時(shí)，會忽略尾流模態(tài)的“慣性”對流場影響[8]。最近，Kim等[9-12]采用數(shù)值模擬方法對圓柱渦激振動開展研究，更深入闡述渦激振動的機(jī)理。本文在前人的基礎(chǔ)上利用開源軟件OpenFOAM采用勻加速，勻減速和勻速三種初始條件對質(zhì)量比為2.6雙自由圓柱渦激振動模擬。研究考慮遲滯效應(yīng)與不考慮遲滯時(shí)，圓柱振動幅值，鎖定區(qū)間，升力系數(shù)，阻力系數(shù)以及尾渦的變化。同時(shí)針對遲滯機(jī)理進(jìn)行討論，并研究了加速度的大小對遲滯現(xiàn)象的影響。

1 計(jì)算模型與模擬工況

1.1 計(jì)算模型建立

本文基于OpenFoam嵌入自編程序進(jìn)行圓柱渦激振動數(shù)值模擬，計(jì)算域描述如圖1所示。

圖1 計(jì)算域幾何示意圖Fig.1 Digram of computational domain

計(jì)算區(qū)域?yàn)樽筮吋吧舷逻吘鄨A心6D，右邊距離圓心20D的矩形流域，采用URANS二維模型，在z向設(shè)置為“empty”邊界條件。網(wǎng)格劃分必須確?！暗屠字Z數(shù)”處內(nèi)的網(wǎng)格質(zhì)量，“低雷諾數(shù)”并不是指整體雷諾數(shù)較低，而是指黏性子層處的雷諾數(shù)較低，為了描述黏性子層內(nèi)的流動引入無量綱數(shù)y+，y+與壁面剪切應(yīng)力有關(guān)。y+的計(jì)算公式由式(1)得出

(1)

式中:τw為壁面剪切力。

網(wǎng)格劃分必須確保第一層網(wǎng)格厚度能使得y+<1。為了保證圓柱面處網(wǎng)格的厚度同時(shí)保證整體網(wǎng)格不過多導(dǎo)致計(jì)算時(shí)間冗長，因此在網(wǎng)格劃分時(shí)，需要將圓柱周圍局部加密，將計(jì)算域劃分為9個(gè)區(qū)域，圓柱周圍4D的范圍內(nèi)采用O型網(wǎng)格。遠(yuǎn)場區(qū)域采用結(jié)構(gòu)網(wǎng)格，其尺寸逐步遞增，如圖2 所示。

采用SST(Shear Stress Transport)k-ω湍流模型和非穩(wěn)態(tài)一階隱式進(jìn)行求解，動量方程的壓力速度耦合采用PIMPLE算法，動量、湍流動能、耗散率項(xiàng)均采用二階迎風(fēng)格式。

(a) (b)圖2 網(wǎng)格計(jì)算區(qū)域劃分示意圖Fig.2 Grid computing region division

1.2 時(shí)間步長控制以及網(wǎng)格無關(guān)性

時(shí)間步長影響著收斂速度和穩(wěn)定性。合理的調(diào)整時(shí)間尺度，能夠在獲得理想收斂速度和穩(wěn)定性的同時(shí)縮短計(jì)算時(shí)間。庫朗數(shù)(Courant Number)指的是時(shí)間步長和物理空間尺度的對應(yīng)關(guān)系，對于湍流而言有著重要影響，庫朗數(shù)由式(2)表示

Cn=u·Δt/Δx

(2)

式中:Cn為庫朗數(shù);u為來流速度;Δt為時(shí)間步長;Δx為網(wǎng)格尺寸。為滿足Courant-Friedrichs-Lewy條件,必須要求庫朗數(shù)足夠小。使用數(shù)值模擬計(jì)算時(shí)候，設(shè)定庫朗數(shù)最大值小于0.50，時(shí)間步長取為Δt=0.000 1。

通過Re=104的圓柱繞流為例驗(yàn)證網(wǎng)格無關(guān)性，采用第一層網(wǎng)格厚度為0.02 mm，近壁面第一層網(wǎng)格的長厚比設(shè)置為8∶1。根據(jù)網(wǎng)格加密與放疏原則，設(shè)置了不同的五種網(wǎng)格方案，不同網(wǎng)格方案的詳細(xì)參數(shù)如表1所示。

表1 不同網(wǎng)格設(shè)計(jì)方案

將阻力與斯特勞哈爾數(shù)的結(jié)果隨網(wǎng)格總體數(shù)量的關(guān)系繪制在圖3和圖4中，由圖3和圖4可知，當(dāng)網(wǎng)格數(shù)達(dá)到以及超過51 502時(shí)，計(jì)算的阻力系數(shù)與斯特勞哈爾數(shù)均達(dá)到了比較好的收斂。而網(wǎng)格數(shù)低于51 502時(shí)，阻力與斯特勞哈爾數(shù)偏差較大，由圖3可以知：網(wǎng)格3已經(jīng)能夠滿足計(jì)算精度的需求。

圖3 阻力均值隨網(wǎng)格數(shù)量變化趨勢Fig.3 The mean drag coefficient with different grid

圖4 斯特勞哈爾數(shù)隨網(wǎng)格數(shù)量變化趨勢Fig.4 The mean drag coefficient with different grid

1.3 模擬工況

選取Javutis等(2004)所做的經(jīng)典的m*=2.6的雙自由度實(shí)驗(yàn)作為參考，圓柱直徑為0.038 1 m，質(zhì)量阻力比m*ξ=0.13，固有頻率為0.4 Hz，流向與橫向頻率比fnx/fny=1,約化速度范圍為2～14，約化速度間隔為0.2。

本文對速度的控制是在OpenFOAM的速度自定義函數(shù)中實(shí)現(xiàn)的，在數(shù)值模擬過程中，三種約化速度的實(shí)現(xiàn)方式如下定義：

(1)勻加速。約化速度從2逐步增加，加速度為每無量綱單元時(shí)間(Ut/D)約化速度增加0.025，增加至目標(biāo)速度后保持該約化速速繼續(xù)模擬30個(gè)周期，同時(shí)上一個(gè)折合流速模擬結(jié)束時(shí)刻的物理場是下一個(gè)折合流速模擬的初始條件。

(2)勻減速。約化速度從14逐步減少，加速度為每無量綱單元時(shí)間(Ut/D)約化速度減小0.025，減小至目標(biāo)速度后保持該約化速速繼續(xù)模擬30個(gè)周期，同時(shí)上一個(gè)折合流速模擬結(jié)束時(shí)刻的物理場是下一個(gè)折合流速模擬的初始條件。

(3)勻速。在數(shù)值模擬開始時(shí)，約化速度直接給定為目標(biāo)約化速度，隨后保持該約化速度繼續(xù)模擬。

本文對速度的控制是在OpenFOAM的速度自定義函數(shù)中實(shí)現(xiàn)的，三種初始條件下圓柱渦激振動各工況參數(shù)如表2所示。

表2 三種初始條件下數(shù)值模擬工況

2 計(jì)算結(jié)果與分析

2.1 渦激振動幅值響應(yīng)

圖3給出了不同初始條件下雙自由度圓柱渦激振動橫向振幅隨約化速度變化曲線，方塊代表Jauvtis等的試驗(yàn)結(jié)果，圓點(diǎn)代表采用勻加速方法，正三角代表采用勻減速方法，倒三角代表采用勻速方法。由圖5可知，三種方式均捕獲到初始分支和下端分支，采用勻加速方式時(shí)，在Ur(約化速度)=6.8時(shí)，Ay/D達(dá)到最大值1.4D。在Jauvtis等的m*=2.6試驗(yàn)中，Jauvtis等發(fā)現(xiàn)在約化速度為8.3，橫向達(dá)到最大振幅比1.5D，他定義該分支為“超長端分支”，因此可知勻加速數(shù)值模擬可以捕獲到“超上端分支”現(xiàn)象，但是其值較為提前并且結(jié)果略小。當(dāng)Ur>6.8時(shí)，振幅下降至0.75D步入下端分支。采用勻減速方式時(shí)，捕獲到上端分支，上端分支與下端分支交界點(diǎn)為Ur=6時(shí)，其最大橫向振幅為0.96D，當(dāng)約化速度大于6時(shí)，振幅比Ay/D陡然下降到0.7D附近，進(jìn)入“下端分支”。采用固定速度方式時(shí)，沒有捕獲到上端分支，初始分支與下端分支交界點(diǎn)為Ur=5，其最大振幅0.85D。在下端分支中，勻加速的振幅略大于試驗(yàn)而勻減速和勻速的振幅略小于實(shí)驗(yàn)值。采用勻加速與勻減速初始條件可以數(shù)值模擬出遲滯現(xiàn)象，如圖5中箭頭所示。然而值得注意的是，Jauvits等(2004)試驗(yàn)的遲滯區(qū)域?yàn)?.8～8.6，而數(shù)值模擬得到的遲滯區(qū)域?yàn)?～6.8。數(shù)值模擬的遲滯區(qū)域相對于試驗(yàn)提前了2左右，但是遲滯區(qū)域的范圍大致都為0.8個(gè)約化速度。

圖5 不同初始條件下圓柱橫向最大振幅Fig.5 Transverse max amplitude in different initial conditions

圖6給出了不同初始條件下圓柱渦激振動流向振幅隨約化速度變化曲線。由圖6可知，采用勻加速方式時(shí)，在Ur=6.8時(shí)，Ax/D達(dá)到最大值0.29D。Jauvtis等的試驗(yàn)中在約化速度為7.3時(shí)達(dá)到中最大流向振幅比為0.32D，數(shù)值模擬結(jié)果較為提前并且結(jié)果偏小。當(dāng)Ur>6.8時(shí)，振幅下降至0.07D步入下端分支。采用勻減速方式時(shí)，在Ur=6時(shí)，Ax/D達(dá)到最大值0.11D，數(shù)值模擬得到的遲滯區(qū)域?yàn)?～6.8。

對比圖5與圖6可知，當(dāng)Ur≤4時(shí)，渦激振動處于初始分支中，圓柱體兩側(cè)漩渦對稱地脫落，產(chǎn)生的升力互相抵消一部分，此時(shí)橫流向運(yùn)動雜亂且運(yùn)動幅值較小，而順流向運(yùn)動規(guī)則且幅值較大，軌跡呈扁平狀。在初始分支內(nèi)，三種速度方式的流向振幅均保持較高的振動幅值。當(dāng)Ur>4時(shí)，渦激振動進(jìn)入鎖定區(qū)域，圓柱形圓柱體兩側(cè)漩渦開始交替脫落，此時(shí)，橫向作用力使橫向運(yùn)動幅值變大，同時(shí)順流向運(yùn)動幅值也增大。進(jìn)入鎖定區(qū)后，采用勻加速與勻減速的流向振幅隨Ur的增大而逐步增大至極值。采用勻速初始條件時(shí)，由于流向振幅進(jìn)入鎖定區(qū)后沒有捕獲到上端分支，因此流向振幅表現(xiàn)為“雙峰”現(xiàn)象，在初始分支Ur=3時(shí)達(dá)到極值0.14D，在Ur=8時(shí)達(dá)到次峰0.075D。

圖6 不同初始條件下圓柱渦激流向振幅Fig.6 In-line max amplitude in different initial conditions

2.2 無因次約化振動頻率響應(yīng)

圖7給出了圓柱渦激振動橫向約化頻率隨約化速度變化曲線，在圖中添加兩條輔助線，實(shí)線為運(yùn)動頻率等于自振動頻率輔助線，虛線為固定圓柱卸渦頻率預(yù)估線。由圖7可知，采用勻加速條件時(shí)：圓柱無無因次鎖定頻率在1.3附近。采用勻減速和勻速條件的鎖定頻率在1.25附近，數(shù)值模擬的鎖定頻率均小于實(shí)驗(yàn)值的1.35。三種初始條件下的鎖定區(qū)域有很大的不同。采用勻加速條件時(shí)，渦激振動鎖定區(qū)范圍為4.2～12.8，其中處于上端分支的為4.2～6.8，處于下端分支的為6.8～12.8；當(dāng)采用勻減速條件時(shí)，渦激振動鎖定區(qū)范圍也為4.2～12.8，但處于上端分支的為4.2～6，處于下端分支的為6～12.8；采用勻速條件時(shí)，渦激振動鎖定范圍為3.8～14，其中處于上端分支的為3.8～5，處于下端分支的為5～14。在鎖定區(qū)間內(nèi)圓柱保持高幅值振動，在鎖定區(qū)外，振動頻率離開實(shí)線向虛線靠攏，表現(xiàn)為固定圓柱泄渦頻率。采用勻加速與勻減速可以得到鎖定區(qū)遲滯區(qū)為6～6.8，如圖7中箭頭所示，其相對于試驗(yàn)遲滯區(qū)7.8～8.6提前約2個(gè)約化速度。

圖7 不同初始條件下圓柱橫向振動約化頻率對比圖Fig.7 Transverse reduced frequency in three initial conditions

采用勻加速初始條件時(shí)，圓柱振動頻率功率函數(shù)密度分布如圖8所示，當(dāng)Ur=3時(shí)，橫向振動頻率等于固有圓柱泄渦頻率，流向振動具有兩個(gè)主要頻率，較大振動頻率對應(yīng)于兩倍橫向振動頻率，較小頻率對應(yīng)單倍橫向頻率，觀察不同約化速度下振動頻率可知，當(dāng)約化速度較低時(shí)(Ur≤4.2)，流向振動出現(xiàn)了兩個(gè)主要頻率，分別為單倍與雙倍橫向頻率。當(dāng)約化速度Ur=6時(shí)，圓柱振動處于鎖定區(qū)內(nèi)，圓柱橫向振動頻率集中在一個(gè)頻率上即圓柱固有頻率，值得注意的是，在鎖定區(qū)域內(nèi)，發(fā)現(xiàn)在除了主頻外，還存在一個(gè)3倍的低頻分量。Jauvtis等(2004)認(rèn)為，這是由于超上端分支產(chǎn)生的“2T”泄渦模型從而導(dǎo)致產(chǎn)生三階升力分量。流向振動頻率為橫向頻率2倍，同時(shí)流向頻率還存在一個(gè)2倍的低頻分量。當(dāng)Ur增大到12.6時(shí)，圓柱處于鎖定區(qū)向解鎖區(qū)跳躍的過渡區(qū)，圓柱橫向振動有兩個(gè)主要頻率，前一個(gè)頻率對應(yīng)固有頻率，后一頻率對應(yīng)于固定圓柱泄渦頻率。流向主要頻率為橫向頻率的2倍，同時(shí)出現(xiàn)1.5倍和2倍的低頻分量。

在鎖定區(qū)域內(nèi)，橫向振動除了主頻外，還存在一個(gè)3倍的低頻分量，說明在鎖定區(qū)域內(nèi)升力存在3倍低頻分量，圖9給出了Ur=6時(shí)渦激升力功率譜密度曲線。由圖9可知：采用勻加速條件下，鎖定區(qū)域內(nèi)渦激升力除主頻外，還存在3倍低頻以及微弱5倍低頻。在鎖定區(qū)域內(nèi)，橫向振動受3倍低頻升力影響從而產(chǎn)生3倍振動分量。

(a)Ur=3 (b)Ur=6 (c)Ur=12.6圖8 各約化速度下橫向振動功率譜密度曲線Fig.8 Transverse vibration frequency PSD in different initial conditions

圖9 Ur=6時(shí)渦激升力功率譜密度曲線Fig.9 Lifting force frequency PSD in increase initial conditions

鎖定區(qū)域升力3倍低頻分量可由泄渦模式所解釋，在鎖定區(qū)域內(nèi)泄渦模型為“2T”模式，如圖8所示。每個(gè)“T”為一個(gè)正負(fù)交替的三重渦，一個(gè)“T”模式可以看為一對半的“2S”泄渦，如圖10圓圈中為一對“2S”泄渦，認(rèn)為“2T”泄渦模式產(chǎn)生了升力中的3倍低頻。

圖10 “T”模式類似于一對半的“2S”模式Fig.10 “T” mode is similar to one and a half of the “2S” mode

2.3 流體力特性響應(yīng)

勻加速，勻減速以及勻速三種不同速度形式所得到的升力均方根和阻力均值隨約化速度變化規(guī)律如圖11所示，對于升力系數(shù)而言，采用勻加速方式時(shí)，在2

圖11 三種速度形式下升力均方根與阻力均值隨約化速度變化圖Fig.11 The RMS lift and mean drag coefficients versus Ur in different initial conditions

對比三種初始條件下的升力系數(shù)與阻力系數(shù)結(jié)果可知，作用在圓柱表面的水動力系數(shù)也存在與圓柱振幅響應(yīng)相類似的遲滯區(qū)域，其區(qū)間范圍為6≤Ur≤6.8。在遲滯區(qū)間內(nèi)，采用勻加速條件的升力系數(shù)較大，其對應(yīng)的橫向振幅也較大。反之，采用勻減速條件的升力系數(shù)較小，其對應(yīng)的橫向振幅也較小。

圖12給出了采用勻加速初始條件時(shí)無量綱位移y/D升力系數(shù)Cl以及阻力系數(shù)Cd隨時(shí)間變化曲線，當(dāng)Ur=3時(shí)，此時(shí)位于“初始分支”，升力系數(shù)與阻力系數(shù)以及橫向位移均表現(xiàn)為規(guī)則的簡弦形式，升力曲線與位移曲線的波峰波谷處在相同位置，表現(xiàn)為典型的同相位角形式。當(dāng)Ur=4時(shí)，圓柱上的升力以及阻力系數(shù)以及橫向位移曲線出現(xiàn)了“拍”現(xiàn)象，此區(qū)域?yàn)槌跏挤种蛏隙朔种н^渡階段，說明當(dāng)圓柱振動從非鎖定區(qū)進(jìn)入鎖定區(qū)時(shí)，存在兩個(gè)相近頻率。當(dāng)Ur=5時(shí)，圓柱振動由初始分支跳躍到上端分支，升力系數(shù)與阻力系數(shù)以及橫向運(yùn)動幅值均有著較大的增加，橫向位移與升力系數(shù)處于同相位。當(dāng)Ur=7時(shí)，圓柱振動從上端分支跳躍至下端分支了，由圖12(d)可知，進(jìn)入下端分之后，升力系數(shù)與阻力系數(shù)以及位移均大幅度下降，升力系數(shù)與橫向振幅曲線由同相位變?yōu)榉聪辔?，出現(xiàn)“相位切換”現(xiàn)象。當(dāng)Ur=9時(shí)，升力系數(shù)山腰處出現(xiàn)“小尖峰”現(xiàn)象，并且升力峰值左邊的小尖峰要高于右側(cè)的小尖峰，這是由于鎖定區(qū)存在三倍頻率分量導(dǎo)致的。當(dāng)Ur=13時(shí)，圓柱脫離鎖定區(qū)，升力系數(shù)與阻力系數(shù)以及橫向振幅均急劇減小，阻力系數(shù)幾乎減小為0，升力與振動處于反相位角的情形。

(a)Ur=3 (b)Ur=4 (c)Ur=5

(d)Ur=7 (e)Ur=9 (f)Ur=13圖12 典型約化速度下升力與阻力隨時(shí)間變化曲線Fig.12 The lift and drag coefficients versus time in different Ur

圖13給出了采用勻減速初始條件時(shí)無量綱位移y/D升力系數(shù)Cl以及阻力系數(shù)Cd隨時(shí)間變化曲線，當(dāng)Ur=13時(shí)，圓柱在鎖定區(qū)外，升力系數(shù)較小，阻力系數(shù)幾乎減小為0，并且升力與振動處于反相位角。當(dāng)Ur=8時(shí)，圓柱下端分支中，升力系數(shù)依舊在山腰處出現(xiàn)“小尖峰”現(xiàn)象，升力與振動處于反相位角。當(dāng)Ur=6時(shí)，圓柱由下端分支跳躍至上端分支，升力驟然增大同時(shí)阻力小幅度增大，升力與振幅由反向變?yōu)橥?，出現(xiàn)“相位切換”現(xiàn)象。當(dāng)Ur=5時(shí)，圓柱處于上端分支，升力系數(shù)與阻力系數(shù)增大，升力與位移表現(xiàn)為同相。Ur=4時(shí)，處于“初始分支”向“上端分支”過渡區(qū)域內(nèi)，升力以及阻力系數(shù)出現(xiàn)了“拍”現(xiàn)象，當(dāng)Ur=3時(shí)，升力系數(shù)與阻力系數(shù)回歸簡諧形式，升力與位移表現(xiàn)為同相。

(a)Ur=13 (b)Ur=8 (c)Ur=6

(d)Ur=5 (e)Ur=4 (f)Ur=3圖13 典型約化速度下升力與阻力以及橫向位移隨時(shí)間Fig.13 The lift and drag coefficients versus time in different Ur

2.4 從運(yùn)動軌跡分析渦激振動

圖14給出采用勻加速方式時(shí)圓柱中心運(yùn)動軌跡，圖中橫坐標(biāo)為流向約化位移，縱坐標(biāo)為橫向約化位移。當(dāng)Ur=2時(shí)，主要發(fā)生流向振動，橫向振動幅值很小。當(dāng)Ur=3時(shí)，圓柱橫向運(yùn)動幅值增大，流向振動頻率變?yōu)闄M向頻率兩倍，運(yùn)動軌跡表現(xiàn)為扁平的“8”字型。當(dāng)Ur=5時(shí)，圓柱進(jìn)入上端分支，此時(shí)橫向與流向的運(yùn)動頻率較為單一，圓柱表現(xiàn)為瘦長的8字型。當(dāng)Ur=6.8時(shí)，圓柱振動軌跡表現(xiàn)出了超上端分支的特征。圓柱橫向振動大幅增大，8字型運(yùn)動軌跡被拉扯為月牙形。當(dāng)Ur=8時(shí)，橫向振動幅值急劇減小至0.7D，圓柱軌跡變?yōu)槭蓍L的8字型。當(dāng)約化速度達(dá)到13.6時(shí)，此時(shí)圓柱振動已經(jīng)脫離鎖定區(qū)，圓柱橫向振幅僅剩0.1D左右。

圖15給出了三種速度方式下各約化速度時(shí)的運(yùn)動軌跡，通過觀察各約化速度下的運(yùn)動軌跡可以清晰地認(rèn)識到低質(zhì)量比雙自由度圓柱渦激振動的振動形式衍化過程。對比三種方式可知，采用勻加速的數(shù)值模擬方式能夠很好的模擬出初始分支，上端分支，超上端分支以及下端分支的運(yùn)動軌跡現(xiàn)象，其橫向振幅呈現(xiàn)遞增趨勢直至Ur=6.8，當(dāng)Ur>6.8時(shí)，其軌跡突變?yōu)槭蓍L8字型。采用勻減速與勻速軌跡差異主要表現(xiàn)在5≤Ur≤6.8之間，勻減速軌跡在Ur>6時(shí)便突變?yōu)槭蓍L8字型。采用勻速方式時(shí)，在Ur>5時(shí)便突變?yōu)槭蓍L8字型。

(a)Ur=2 (b)Ur=3 (d)Ur=5

(f)Ur=6.8 (g)Ur=8 (h)Ur=14圖14 勻加速各約化速度下的軌跡圖Fig.14 The VIV trajectory versus time in different Ur

圖15 不同初始條件下圓柱渦激振動軌跡圖Fig.15 The VIV trajectory versus time in different initial conditions

2.5 從渦脫形式分析遲滯現(xiàn)象

圖16給出了勻加速初始條件下圓柱泄渦等值云圖(左為Jauvtis等試驗(yàn)DPIV結(jié)果，右為CFD結(jié)果)。Ur=2時(shí)，處于初始分支中，每半個(gè)周期便會在圓柱一邊卸載一個(gè)漩渦，該泄渦模式為“2S”模式。Ur=5，處于上端分支中，每半個(gè)周期便會在圓柱一邊卸載一對漩渦，同時(shí)一個(gè)漩渦的渦量強(qiáng)度明顯大于另一個(gè)漩渦，該泄渦模式為“2P”模式。Ur=6.8，處于超上端分支中，每半個(gè)周期便會在圓柱一邊卸載三個(gè)漩渦，該泄渦模式為“2T”模式，當(dāng)圓柱橫向得到超大振幅時(shí)，便會體現(xiàn)出該種泄渦模式。Ur=8時(shí)，每半個(gè)周期便會在圓柱一邊卸載一對漩渦，這兩個(gè)漩渦強(qiáng)度相當(dāng)，該泄渦模式為“2P”模式。采用勻加速初始條件成功模擬出“2T”模式，這得益于考慮了渦激振動“遲滯效應(yīng)”。

圖17給出了勻減速初始條件下圓柱泄渦等值云圖，Ur=6.8時(shí)，泄渦模式表現(xiàn)為“2P”模式，對比圖12(c)可知，當(dāng)Ur=6.8時(shí)，采用勻加速與勻減速兩種初始條件下，圓柱尾渦模式不同，具體為：采用勻加速時(shí)，其尾渦為“2T”模式，與其較小Ur下的模式相同，而采用勻減速時(shí)，其尾渦為“2P”模式，與其較大Ur下的模式相同。勻減速Ur=6時(shí)，此時(shí)流體力超過尾渦維持“2P”模式的最大承受能力，其泄渦模式變?yōu)椤?T”模式，該泄渦模式伴隨著較大的升力，從而導(dǎo)致振幅的突然增大。不同的尾渦模式下升力大小不同，從而導(dǎo)致橫向振幅不同以至于產(chǎn)生了“遲滯現(xiàn)象”。圖18給出勻速初始條件下圓柱泄渦等值云圖，Ur=5時(shí)，泄渦模式表現(xiàn)為“2P”模式，Ur=6.8時(shí)，泄渦模式依舊為“2P”模式。采用勻速初始條件時(shí)，圓柱泄渦一直表現(xiàn)為“2P”模式，該模式難以提供較大振幅，因此勻速初始條件未能捕獲上端分支。對比三種初始條件下的尾渦模式可知，圓柱的泄渦模式體現(xiàn)出“慣性”性質(zhì)，尾渦的傾泄模式總有著維持原傾泄模式的趨勢，當(dāng)約化速度提供的作用力大于圓柱維持原泄渦模式的最大承受能力后，圓柱泄渦模式將出現(xiàn)跳躍，不同初始條件對尾渦的影響進(jìn)一步體現(xiàn)在振幅的差異上，同Ur下的不同泄渦模式導(dǎo)致產(chǎn)生“遲滯現(xiàn)象”。

(i)DPIV結(jié)果 (ii)數(shù)值模擬結(jié)果 (i)DPIV結(jié)果 (ii)數(shù)值模擬結(jié)果(a)Ur=3渦量等值云圖 (b)Ur=5渦量等值云圖

(i)DPIV結(jié)果 (ii)數(shù)值模擬結(jié)果 (i)DPIV結(jié)果 (ii)數(shù)值模擬結(jié)果(c)Ur=6.8渦量等值云圖 (d)Ur=8渦量等值云圖圖16 勻加速圓柱泄渦等值云圖對比圖Fig.16 The vorticity magnitude contours versus Ur in increasing condition

(a)Ur=6.8渦量等值云圖(b)Ur=6渦量等值云圖(a)Ur=5渦量等值云圖(b)Ur=6.8渦量等值云圖圖17 勻減速數(shù)值模擬圓柱泄渦等值云圖 圖18 勻速數(shù)值模擬圓柱泄渦等值云圖Fig.17 The vorticity magnitude contours versus Ur in Fig.18 The vorticity magnitude contours versus Ur in decreasing condition using CFD constant condition using CFD

2.6 加速度a對數(shù)值結(jié)果影響

由上節(jié)可知，不同初始條件下的渦激振動體現(xiàn)出很大的差異性，初始條件對圓柱泄渦模式有很大影響，當(dāng)約化速度提供的作用力大于圓柱維持原泄渦模式的最大承受能力后，圓柱泄渦模式將出現(xiàn)跳躍，同時(shí)圓柱振動振幅將出現(xiàn)跳躍，以至于產(chǎn)生“遲滯現(xiàn)象”。增大速度時(shí)，圓柱泄渦將傾向于速度較小的泄渦形式。減小速度時(shí)，圓柱泄渦將傾向于速度較大的泄渦形式。由于流體具有慣性效應(yīng)，因此采用勻加速和勻減速進(jìn)行數(shù)值模擬研究是十分有必要的，本節(jié)將探討采用勻加速初始條件下，加速度對結(jié)果的影響。分別采用0.012 5，0.025，0.05以及0.1的約化加速度進(jìn)行數(shù)值模擬，約化速度范圍為3～10，約化速度間隔為0.5，同時(shí)Ur=6.8為臨界點(diǎn)，也加以研究。

不同加速度下圓柱橫向振幅變化規(guī)律如圖19所示。采用四種加速度下，均捕獲到初始分支上端分支與下端分支。采用約化加速度為0.012 5，0.025的方式捕獲到了超上端分支。在3≤Ur≤5.5時(shí)，四種加速度方式得到的振幅基本一致，在該分支中，加速度對振動形態(tài)影響不大。在5.5

圖19 不同加速度下圓柱渦激橫向振幅Fig.19 Transverse max amplitude with different acceleration

圖20給出了不同加速度下升力均方根值與阻力均值變化曲線。在3≤Ur≤5.5時(shí)，四種加速度得到的升力基本一致。在5.5

(a)

(b)圖20 四種加速度下升力均方根與阻力均值隨約化速度關(guān)系Fig.20 The RMS lift coefficient and mean drag coefficient with different accelerations

圖21給出了不同加速度下遲滯區(qū)間變化曲線。采用約化加速度為0.012 5，0.025的方式捕獲到了遲滯區(qū)間，其范圍為6≤Ur≤6.8。采用約化加速度為0.05，0.1的方式?jīng)]有捕獲到超上端分支，其遲滯區(qū)間為5.5≤Ur≤6。由此可見，加速度大小對遲滯區(qū)間有較大影響，約化加速度在0.025～0.05存在一個(gè)臨界值，小于臨界值與大于臨界值的渦激振動遲滯區(qū)域不同。但同在大于臨界值或同在小于臨界值時(shí)，渦激振動得到遲滯區(qū)域相同。

圖21 不同加速度下渦激振動遲滯區(qū)域Fig.21 Hysteresis region of VIV under different acceleration

圖22給出Ur從6.5加速到6.8時(shí)升力阻力以及位移變化曲線，其采用約化加速度分別為0.025和0.05?？梢钥吹?，約化加速度為0.025時(shí)，升力系數(shù)與阻力系數(shù)隨著時(shí)間逐漸增加，而采用約化加速度為0.05時(shí)，升力系數(shù)與阻力系數(shù)隨著時(shí)間增加而陡然下降，同時(shí)對比橫向位移與升力系數(shù)可知，此時(shí)發(fā)生了相位切換。這說明采用較大加速度將會提前導(dǎo)致“相位切換”，從而使得振幅減小。

(a)a*=0.025

(b)a*=0.05圖22 Ur=6.8時(shí)下升力與阻力隨時(shí)間變化曲線Fig.22 The lift coefficient and drag coefficient versus time

3 結(jié) 論

本文基于PimpleDYMFoam求解器結(jié)合SSTk-ω湍流模型，進(jìn)行了雙自由度低質(zhì)量數(shù)值模擬研究，模擬采用了三種初始條件：勻加速，勻減速以及勻速條件。數(shù)值模擬捕獲到遲滯現(xiàn)象并且三種條件下圓柱振動響應(yīng)有顯著區(qū)別，數(shù)值模擬對比得出結(jié)果如下：

(1)在低質(zhì)量比圓柱渦激振動試驗(yàn)中，采用加速與減速兩種試驗(yàn)條件時(shí)，上端分支與下端分支存在遲滯。在“數(shù)值水池”中采用勻加速與勻減速方式也可捕獲遲滯現(xiàn)象。同時(shí)發(fā)現(xiàn)，采用勻加速與勻減速條件能捕獲到超上端分支，而采用勻速條件則缺失上端分支。在Ur=6.8時(shí)，勻加速方式得到最大振幅1.4D，而采用不考慮遲滯效應(yīng)的勻速條件時(shí)，圓柱最大振幅僅為0.85D。由此可知，遲滯效應(yīng)的存在對振幅的大小以及突變具有顯著影響。

(2)遲滯效應(yīng)的機(jī)理可以通過泄渦模式來解釋，圓柱泄渦模式體現(xiàn)出“慣性”性質(zhì)，尾渦的傾泄模式總有著維持原傾泄模型的趨勢，即增大速度時(shí)，圓柱泄渦將傾向于速度較小的泄渦形式。減小速度時(shí)，圓柱泄渦將傾向于速度較大的泄渦形式。當(dāng)約化速度提供的作用力大于圓柱維持原泄渦模式的最大承受能力后，圓柱泄渦模式將出現(xiàn)跳躍。在遲滯區(qū)間內(nèi)，同Ur下采用不同初始條件時(shí)的泄渦模式是不同的，同Ur下不同的泄渦模式導(dǎo)致了“遲滯現(xiàn)象”。

(3)采用勻加速勻減速以及勻速條件均捕獲到“鎖定”現(xiàn)象，但是鎖定區(qū)間不同。上端分支跳躍至下端分支的約化速度大小與初始條件有關(guān)。同時(shí)發(fā)現(xiàn)，當(dāng)Ur≤4.2時(shí)，流向振動有兩個(gè)主要頻率，分別為單倍與雙倍橫向頻率。在鎖定區(qū)域內(nèi)，橫向振動除主頻外，還存在一個(gè)3倍的低頻分量。

(4)作用在圓柱體上的升力與阻力也出現(xiàn)了遲滯現(xiàn)象，其遲滯區(qū)域?yàn)?≤Ur≤6.8，并在約化速度為6和6.8時(shí)，通過升力與位移關(guān)系捕獲了相位切換現(xiàn)象。升力與橫向位移剛開始為同向，隨著流速增加，泄渦模式發(fā)生突變導(dǎo)致圓柱受到的升力發(fā)生突變，從而導(dǎo)致升力與橫向位移變?yōu)榉聪?，認(rèn)為是泄渦模式的突變導(dǎo)致了相位切換。

(5)勻加速條件出現(xiàn)了“月牙形”運(yùn)動軌跡，而勻減速與勻速只捕獲了經(jīng)典“8”型運(yùn)動。同時(shí)發(fā)現(xiàn)，在遲滯區(qū)間內(nèi)，相同約化速度下的運(yùn)動軌跡形式是不同的，圓柱軌跡也與初始條件有一定關(guān)系。

(6)在遲滯區(qū)域內(nèi)，在約化加速度0.025與0.05間存在臨界值，大于該臨界值，升力與位移會提前發(fā)生相位切換，小于該臨界值，相位切換將會延后。數(shù)值模擬中，控制約化加速度小于等于0.025是其捕獲超上端分支的一個(gè)必要條件。

[ 1 ] 白治寧, 肖龍飛, 程正順,等. 深吃水半潛式平臺渦激運(yùn)動響應(yīng)模型實(shí)驗(yàn)研究[J]. 船舶力學(xué), 2014(4):377-384.

BAI Zhining, XIAO Longfei, CHENG Zhengshun, et al. Experimental study on vortex induced motion response of a deep draft semi-submersible platform[J]. Journal of Ship Mechanics， 2014(4): 377-384.

[ 2 ] 曹淑剛, 黃維平, 顧恩凱. 考慮流固耦合的彈性圓柱體渦激振動研究[J]. 振動與沖擊, 2015, 34(1): 58-62.

CAO Shugang, HUANG Weiping, GU Enkai. Vortex-induced vibration of an elastic cylinder considering fluid-structure interaction[J]. Journal of Vibration & Shock, 2015, 34(1): 58-62.

[ 3 ] 黃維平, 劉娟, 唐世振. 考慮流固耦合的大柔性圓柱體渦激振動非線性時(shí)域模型[J]. 振動與沖擊, 2012, 31(9): 140-143.

HUANG Weiping, LIU Juan, TANG Shizhen. Nonlinear model of vortex induced vibration of flexible cylinder in consideration of fluid-structure interaction[J]. Journal of Vibration & Shock, 2012, 31(9): 140-143.

[ 4 ] 范杰利, 黃維平. 細(xì)長立管兩向自由度渦激振動數(shù)值研究[J]. 振動與沖擊, 2012, 31(24): 65-68.

FAN Jieli, HUANG Weiping. Numerical simulation of 2-DOF vortex-induced vibration of a long riser[J]. Journal of Vibration & Shock, 2012, 31(24): 65-68.

[ 5 ] JAUVTIS N，WILLIAMSON C H K. The effect of two degrees of freedom on vortex-induced vibration at low mass and damping[J]. Journal of Fluid Mechanics, 2004, 509(509): 23-62.

[ 6 ] 谷家揚(yáng), 楊建民, 肖龍飛. 兩種典型立柱截面渦激運(yùn)動的分析研究[J]. 船舶力學(xué), 2014, 18(10): 1184-1194.

GU Jiayang, YANG Jianmin, XIAO Longfei. Study on vortex induced motion of two typical different cross-section columns[J]. Journal of Ship Mechanics, 2014, 18(10): 1184-1194.

[ 7 ] GUILMINEAU E, QUEUTEY P. Numerical simulation of vortex-induced vibration of a circular cylinder with low mass-damping in a turbulent flow[J]. Journal of Fluids and Structures, 2004, 19(4): 449-466.

[ 8 ] 陳威霖, 及春寧, 徐萬海. 并列雙圓柱流致振動的不對稱振動和對稱性遲滯研究[J]. 力學(xué)學(xué)報(bào), 2015, 47(5): 731-739.

CHEN Weilin, JI Chunning, XU Wanghai. Numerical investigation on the asymmetric vibration and symmetry hysteresis of flow-induced vibration of two side-by-side cylinders[J]. Chinese Journal of Theoretical & Applied Mechanics, 2015, 47(5): 731-739.

[ 9 ] KIM S, WILSON P A, CHEN Z M. Numerical simulation of force and wake mode of an oscillating cylinder[J]. Journal of Fluids & Structures, 2014, 44(7): 216-225.

[10] KHALAK A, WILLIAMSON C H K. Motions, forces and mode transitions in vortex-induced vibrations at low mass-damping[J]. Journal of Fluids & Structures, 1999, 13(7/8): 813-851.

[11] MINGUEZ M，LUPPI A，BERGER A. Slender buoy FSHR vortex induced rotations[C]// Proceedings of the ASME 31th International Conference on Ocean, Offshore and Arctic Engineering. Rio de Janeiro: OMAE, 2012.

[12] CHEN W, JI C, XU W, et al. Response and wake patterns of two side-by-side elastically supported circular cylinders in uniform laminar cross-flow[J]. Journal of Fluids & Structures, 2015, 55: 218-236.