單自由度圓柱渦激振動特性及其影響因素分析?

袁方正， 袁德奎， 曾 攀， 劉長根

(天津大學(xué)機械工程學(xué)院力學(xué)系，天津 300354)

在來流作用下深水立管發(fā)生渦激振動 (Vortex induced vibration, VIV) 是一個在海洋工程中普遍的現(xiàn)象，同時也是引起立管產(chǎn)生疲勞損壞的主要原因之一。水流流經(jīng)鈍尾體如圓柱類結(jié)構(gòu)時，會在尾流區(qū)發(fā)生流動分離并產(chǎn)生周期性泄放的漩渦，尾渦的周期性泄放會導(dǎo)致圓柱受到周期性的荷載，從而誘發(fā)渦激振動。渦激振動的多自由度、強非線性的特點決定了它是一個復(fù)雜的流固耦合問題。當(dāng)立管結(jié)構(gòu)發(fā)生渦激振動且振動頻率接近系統(tǒng)自身的固有頻率時，會發(fā)生共振，促使結(jié)構(gòu)物產(chǎn)生大幅度的振動，從而對海洋結(jié)構(gòu)物產(chǎn)生破壞。隨著人們對海洋資源的開發(fā)，對圓柱渦激振動現(xiàn)象的研究也顯得更為迫切。

實驗是研究渦激振動問題的主要的方法之一，大多數(shù)對渦激振動特性和其產(chǎn)生機理的認(rèn)識都是來源于實驗。Feng[1]在風(fēng)洞中展開了高質(zhì)量比低阻尼比的單自由度渦激振動實驗，該實驗發(fā)現(xiàn)振動存在兩個響應(yīng)分支，并觀測到分支在切換時存在遲滯(Hysteretic loop)現(xiàn)象。Williamson團(tuán)隊在水中對圓柱的渦激振動進(jìn)行了大量實驗研究[2-5]，發(fā)現(xiàn)圓柱進(jìn)行渦激振動的鎖定區(qū)間范圍與質(zhì)量比m*有關(guān)?，F(xiàn)有的實驗不僅增加了學(xué)者對圓柱渦激振動特征和機理的認(rèn)知，也為檢驗和發(fā)展相應(yīng)的數(shù)值模擬方法奠定了基礎(chǔ)。計算機性能的提高使CFD(計算流體力學(xué))方法成為研究渦激振動的主要方法之一。Sanchis等[6]認(rèn)為質(zhì)量比是影響圓柱體振動特性的重要因素。陳正壽等[7]討論了質(zhì)量比不同的圓柱在相同約化速度下，由于單雙自由度的不同所造成的振動特性的差別，簡要分析了順流向振動對橫流向振動的影響程度。曹淑剛等[8]探討了不同質(zhì)量比對圓柱渦激振動的順流向頻率的影響。此前的研究工作對有關(guān)質(zhì)量阻尼比m*ζ與阻尼比ζ對振動特性影響的研究較少，也沒有從理論，實驗，數(shù)值模擬3個方面綜合考慮這一問題。

本文在已有研究的基礎(chǔ)上，對流體荷載進(jìn)行線性分解推導(dǎo)得到了鎖定區(qū)間柱體響應(yīng)的頻率比、振幅與流體荷載和系統(tǒng)物性參數(shù)之間的定性關(guān)系。根據(jù)推導(dǎo)得出的結(jié)果，結(jié)合相關(guān)的經(jīng)典實驗，對圓柱渦激振動特性的影響因素(質(zhì)量比，阻尼比，質(zhì)量阻尼比)開展了系統(tǒng)的數(shù)值實驗，對比了不同因素下圓柱鎖定區(qū)間、振幅、升阻力系數(shù)等關(guān)鍵參數(shù)和特征，進(jìn)一步分析了不同質(zhì)量比，阻尼比，質(zhì)量阻尼比對圓柱渦激振動特性的影響。

1 線性化分析

基于單自由度剛性圓柱振動的軸向?qū)ΨQ性，可將模型簡化為圖1所示的二維單自由度振動系統(tǒng)。

圖1 單自由度橫向振動的物理模型

彈性支承的二維圓柱單自由度振動系統(tǒng)實際上是一種彈簧-質(zhì)量-阻尼系統(tǒng)，其振動方程為：

(1)

式中：m為單位長度的圓柱質(zhì)量，y為圓柱偏離平衡位置的位移，k為彈簧剛度系數(shù)，c為振動系統(tǒng)的阻尼系數(shù)，F(xiàn)y(t)為圓柱受到的升力。

在鎖定區(qū)，將圓柱的振動近似認(rèn)為簡諧振動，且泄渦頻率完全被固定在圓柱振動頻率上。因此假設(shè)圓柱橫向位移y的表達(dá)式為：

y=Asin(ωext)。

(2)

式中：A為圓柱的振幅；ωex為圓柱的振動圓頻率。圓柱所受升力無量綱化后的升力系數(shù)為：

CL=Clsin(ωext+φ)。

(3)

式中：Cl為升力系數(shù)的幅值；φ表示升力系數(shù)與位移時間歷程的相位角。根據(jù)Parkinson[9]的方法：將升力系數(shù)進(jìn)行分解，如下：

CL=Clvcos(ωext)+Cla[-sin(ωext)]。

(4)

由式(4)和式(5)可得：

Clv=Clsin(φ),Cla=-Clcos(φ)，

(5)

φ=arctan(-Clv/Cla)。

(6)

式中：Clv為橫向流體荷載的阻尼分量，體現(xiàn)系統(tǒng)能量的積累與耗散；Cla為橫向流體荷載的慣性分量，與附加質(zhì)量有關(guān)。

下面對該式(1)進(jìn)行無量綱化：

式中：CL為無量綱升力系數(shù)；Y為無量綱位移，U為來流速度；D為圓柱直徑；U*為約化速度；fex為圓柱實際振動頻率；fn為圓柱固有頻率；ωn為固有圓頻率；ρ為流體密度。將振動方程進(jìn)行無量綱化，可得：

(7)

再將式(2)與式(4)無量綱結(jié)果代入式(2)中，可得：

(8)

(9)

式(8)和式(9)反映了頻率比、振幅與流體荷載和系統(tǒng)物性參數(shù)之間的關(guān)系，可以定性反映出鎖定區(qū)間柱體的振動特性。但是引入了相關(guān)假設(shè)，因此流體荷載系數(shù)分解得到的上述兩式只能在當(dāng)流體荷載頻率完全被圓柱振動頻率所捕獲時才有定性參考作用。由于分析需要，暫且假定當(dāng)圓柱處于鎖定區(qū)間內(nèi)時，流體荷載頻率(泄渦頻率)與振動頻率及圓柱固有頻率都相等(這在高質(zhì)量比是近似成立的，但在低質(zhì)量比情況下會有差異)，并且圓柱的振動可以由正弦函數(shù)來表示。則由式(8) 可知，在鎖定區(qū)間組合因子m*ζ越大，可能導(dǎo)致A*越小；由式 (9) 可見，在鎖定區(qū)間假定頻率比固定，則m*越大可能使得U*越大，進(jìn)入鎖定的約化速度越大即鎖定段延后，鎖定區(qū)間變窄。

2 實驗結(jié)果分析

美國康奈爾大學(xué)的Williamson教授等沿用了Feng的經(jīng)典實驗設(shè)計，對低質(zhì)量比低阻尼比圓柱的渦激振動進(jìn)行了大量實驗研究。根據(jù)上一節(jié)的線性化分析得出的結(jié)論將Williamson系列實驗的結(jié)果進(jìn)行了分類整理(見圖2)。

圖2 Williamson實驗結(jié)果

從圖中可以看出：圓柱振動系統(tǒng)的組合因子m*ζ和質(zhì)量比m*的不同會導(dǎo)致圓柱振動特性出現(xiàn)較大差異。具體表現(xiàn)為組合因子m*ζ決定了圓柱的振幅，在質(zhì)量比一定的情況下，其越小則振動的幅值越大；而在組合因子一定的情況下，圓柱的鎖定區(qū)間范圍主要受質(zhì)量比m*的影響，質(zhì)量比越小則鎖定區(qū)間的范圍越大。這表明線性化分析的結(jié)果與實驗結(jié)果相符，都反映出渦激振動特性受到m*ζ和m*的影響較大。下文將對不同的質(zhì)量比和質(zhì)量阻尼比的圓柱進(jìn)行系統(tǒng)的數(shù)值模擬，以期從數(shù)值模擬的角度上對影響渦激振動特性的因素進(jìn)行定量的分析。

3 數(shù)學(xué)模型及求解方法

3.1 流動控制方程

將問題簡化為不可壓縮流體在與圓柱軸線垂直的平面內(nèi)的二維流動，可寫出如下的雷諾平均納維斯托克斯方程(RANS)：

(10)

3.2 湍流模型

根據(jù)文獻(xiàn)[10]所得出的結(jié)論，采用SSTk-ω湍流模型對RANS方程進(jìn)行封閉。SST(Shear Stress Transport)k-ω模型通過添加一個混合函數(shù)使得模型在近壁面處為k-ω模型，而在遠(yuǎn)離壁面的充分湍動區(qū)則為高雷諾數(shù)k-ε模型，是對k-ω模型的改進(jìn)，輸運方程如下：

(11)

式中：k為湍動能；ω為比耗散率；和分別為k、ω方程的湍流普朗特數(shù)；μt為湍流黏性系數(shù)；F1為混合函數(shù)。模型對湍流黏度的定義考慮了湍流剪切力的輸運，湍流模型的細(xì)節(jié)詳見文獻(xiàn)[11]。該模型考慮了湍流剪切力效應(yīng)，在模擬逆壓梯度、邊界層分離、近壁區(qū)的繞流旋流等有較大優(yōu)勢。

3.3 模型設(shè)置及求解方法

計算網(wǎng)格劃分見圖3。整個計算域長40D、寬20D，圓柱圓心距離計算域上邊界和右邊界離分別為10D和30D。采用混合網(wǎng)格畫法，在圓柱外設(shè)置半徑為4D的圓形區(qū)域，此區(qū)域劃分為結(jié)構(gòu)化網(wǎng)格，使其隨著圓柱一起運動，從而保證了邊界層網(wǎng)格的穩(wěn)定，且適當(dāng)減小由于網(wǎng)格運動產(chǎn)生的網(wǎng)格畸變率導(dǎo)致的數(shù)值不穩(wěn)定，并且對圓柱周邊網(wǎng)格進(jìn)行加密以保證計算的空間分辨率。外部區(qū)域為三角形非結(jié)構(gòu)網(wǎng)格，以保證使用動網(wǎng)格彈性光順技術(shù)時的靈活性。

圖3 計算域網(wǎng)格

利用有限體積法求解渦激振動的流場：時間離散采用二階隱式格式，用二階迎風(fēng)格式對對流項進(jìn)行離散，同時用SIMPLEC算法求解壓力速度耦合方程。用四階Runge-Kutta法求解結(jié)構(gòu)振動方程。每一個時間步內(nèi)的具體求解過程為：(1)求解RANS方程獲得流場；(2)對柱面的應(yīng)力積分以獲得升力；(3)將無量綱化的升力代入結(jié)構(gòu)振動方程中，求解振動方程獲得圓柱的位移和速度；(4)將圓柱的運動信息傳遞給求解器，再用動網(wǎng)格功能實現(xiàn)網(wǎng)格更新。如此反復(fù)迭代，直至達(dá)到設(shè)定的收斂精度，實現(xiàn)流固耦合計算。其中后三步都是基于C語言環(huán)境編寫代碼實現(xiàn)。

3.4 網(wǎng)格收斂檢驗

數(shù)值計算精度和效率很大程度上依賴于網(wǎng)格的質(zhì)量。為兼顧分辨率和計算效率，在開展正式的模擬之前對網(wǎng)格收斂性進(jìn)行檢驗是必需的。本文針對所要開展的模擬工作按以下步驟系統(tǒng)檢驗了網(wǎng)格收斂性。

首先分析了邊界網(wǎng)格(主要是第一層的網(wǎng)格無量綱高度y+和圓周網(wǎng)格數(shù))對計算精度和效率的的影響，并從中選擇最接近實驗結(jié)果的壁面邊界處理方法。參數(shù)設(shè)置參考Williamson[2]的實驗：圓柱直徑D=0.038 1 m，阻尼比ζ=0.005 4，質(zhì)量比m*=2.4，約化速度為9，雷諾數(shù)處于亞臨界范圍。

根據(jù)經(jīng)驗，y+是1的量級，y+要在模擬結(jié)束后才可獲得，因此具體使用時可先采用經(jīng)驗公式[12]：

y+=0.172Re0.9(y1/D)。

(12)

對網(wǎng)格第一層高度進(jìn)行估算，得到結(jié)果后再不斷調(diào)整網(wǎng)格直到達(dá)到符合條件的y+值。本文設(shè)置了兩組工況分別研究了y+、圓周網(wǎng)格數(shù)不同引起結(jié)果的差異，計算結(jié)果見表1。

表1 不同工況下的結(jié)果及對比

A組工況對比了y+的不同，可知y+為3時結(jié)果較理想；B組工況比較了不同的圓周網(wǎng)格數(shù)的影響，主要考察計算效率，可知當(dāng)圓周網(wǎng)格數(shù)為120時，網(wǎng)格滿足收斂要求。綜合來看，采用B2工況進(jìn)行后續(xù)計算可以滿足網(wǎng)格收斂性要求。

4 計算結(jié)果及分析

4.1 不同質(zhì)量比、組合因子對渦激振動特性的影響

為了對渦激振動的影響因素進(jìn)行分析，首先研究了固定阻尼比ζ=0.000 54的情況下，不同質(zhì)量比圓柱渦激振動特性的差異。對比了高低兩種質(zhì)量比工況：m*=2.4，ζ=0.000 54 (m*ζ=0.001 3) 和m*=24，ζ=0.000 54 (m*ζ=0.013)。得出了不同質(zhì)量比圓柱的振幅隨約化速度的變化曲線，結(jié)果見圖4。

圖4 ζ=0.000 54時不同質(zhì)量比圓柱的振幅響應(yīng)

從圖4可以看出當(dāng)固定ζ=0.000 54改變質(zhì)量比時，圓柱振動特性發(fā)生了較大變化。當(dāng)質(zhì)量比減小時，圓柱在鎖定區(qū)的最大振幅略微上升，整體振幅也有較大增長；從鎖定段來看，隨著質(zhì)量比減小，圓柱的鎖定區(qū)間有了大范圍的擴展，大質(zhì)量比情況對應(yīng)鎖定段約為U*∈[4.5,5]，而小質(zhì)量比情況約為U*∈[4,10]，鎖定區(qū)間長度相差了一個數(shù)量級。

圖5給出了不同約化速度下2種質(zhì)量比圓柱振動的力系數(shù)和位移的時間歷程曲線?？梢钥闯觯瑢τ诓煌馁|(zhì)量比，都存在差拍和位移升力相位角180°轉(zhuǎn)變等現(xiàn)象。對于m*=24的圓柱，當(dāng)U*=3時位于初始分支，振幅很?。划?dāng)U*=4，可以觀察到力系數(shù)和位移都出現(xiàn)了明顯的多頻振動的拍現(xiàn)象，表明振動進(jìn)入過渡區(qū)，振幅有所增大；而后振幅再次增大，差拍現(xiàn)象消失，曲線重新變?yōu)榉€(wěn)定簡諧振動狀態(tài)，表明圓柱進(jìn)入鎖定區(qū)；最后當(dāng)U*=6時，位移大幅減小，圓柱進(jìn)入解鎖區(qū)。對于m*=2.4的圓柱，當(dāng)U*=3時進(jìn)入過渡區(qū)；當(dāng)U*=12時，圓柱已經(jīng)進(jìn)入解鎖區(qū)。

圖5 不同約化速度下Cd、Cl、y/D的時間曲線

以上分析表明：在固定ζ=0.000 54的情況下改變質(zhì)量比m*，質(zhì)量比m*和組合因子m*ζ均發(fā)生了變化，由此導(dǎo)致鎖定區(qū)的振幅和鎖定范圍都發(fā)生變化。

為了進(jìn)一步驗證理論分析與實驗給出的結(jié)論，接下來將分別對m*和m*ζ的影響進(jìn)行討論。

4.2 質(zhì)量比對渦激振動特性的影響

在固定組合因子m*ζ=0.013的情況下，選取了m*=1.2、2、2.4、6、24五種質(zhì)量比，通過數(shù)值模擬來研究不同質(zhì)量比對渦激振動特性產(chǎn)生的影響。

圖6給出了質(zhì)量比的變化所導(dǎo)致的渦激振動特性的差異。從圖中可見，隨著質(zhì)量比逐漸減小，圓柱的鎖定范圍有逐漸擴大的趨勢，但在鎖定區(qū)的最大振幅基本相同。此外，可以觀察到鎖定區(qū)間的擴大是“兩頭延伸”的，即低質(zhì)量比比高質(zhì)量比更早進(jìn)入鎖定且退出鎖定的約化速度延遲。從鎖定范圍擴大的位置來看，擴大的鎖定區(qū)主要發(fā)生在高流速一端。當(dāng)質(zhì)量比減小到1.2時，鎖定區(qū)間覆蓋了整個約化速度范圍，并可以觀察到其振幅變化趨勢是逐漸減小而非快速減小到非鎖定段。在圖6同時也表明了振動圓柱的最大振幅由組合因子m*ζ決定，當(dāng)組合因子相等時，圓柱處于鎖定狀態(tài)的最大振動幅值基本相等，鎖定區(qū)間內(nèi)的振幅也基本相等。

圖6 時不同質(zhì)量比圓柱的振幅響應(yīng)

在線性化分析時，由式 (9) 得出的m*越小，鎖定區(qū)間變大這一結(jié)論只在高質(zhì)量比的條件下(假定此時泄渦頻率完全固定在圓柱固有頻率上)成立，而從圖中可以看出，質(zhì)量比較低時，鎖定區(qū)間的范圍仍會隨著質(zhì)量比的減小逐漸變大，這表明線性化分析得出的鎖定區(qū)間范圍的大小與質(zhì)量比成負(fù)相關(guān)關(guān)系的結(jié)論可推廣至低質(zhì)量比的情況。

4.3 組合因子對渦激振動特性的影響

最后，為了研究組合因子的影響，對固定質(zhì)量比m*=2.4并采用不同的阻尼比即不同組合因子的情況進(jìn)行了渦激振動的數(shù)值模擬。

圖7給出了當(dāng)m*=2.4時，圓柱在3種不同組合因子下的渦激振動響應(yīng)。從中可以看出，在固定m*=2.4的情況下，3種組合因子得到的響應(yīng)幅值隨約化速度的變化趨勢是一致的。而且可以看出高低振幅切換的位置也相同，對應(yīng)鎖定區(qū)間的范圍一致。從振幅的角度來看，在鎖定區(qū)間內(nèi)，隨著組合因子的減小，整個區(qū)間的振幅呈現(xiàn)增大趨勢，但是增大的幅度十分有限。脫離鎖定后，3種組合因子得到的解鎖區(qū)振幅響應(yīng)差異不大。

圖7 m*=2.4時不同組合因子圓柱的振幅響應(yīng)

5 結(jié)論

本文在前人實驗結(jié)果的基礎(chǔ)上，采用線性化分析和數(shù)值模擬相結(jié)合的方法，對影響渦激振動特性的參數(shù)進(jìn)行了分析，結(jié)果表明：

(1)鎖定區(qū)間柱體響應(yīng)的頻率比、振幅與流體荷載和系統(tǒng)物性參數(shù)之間具有確定的關(guān)系：圓柱的振幅與質(zhì)量阻尼比成負(fù)相關(guān)關(guān)系，鎖定區(qū)間范圍的大小與質(zhì)量比成負(fù)相關(guān)關(guān)系。

(2)當(dāng)圓柱處在鎖定狀態(tài)時，其振幅大小由組合因子m*ζ決定。組合因子越小，圓柱的振幅越大，反之則越小，但是組合因子對振幅的影響幅度有限；柱體鎖定區(qū)間范圍的大小由質(zhì)量比m*決定，質(zhì)量比越小，鎖定范圍越大，反之則越小，質(zhì)量比對鎖定區(qū)間范圍的大小有顯著影響，影響的主要范圍是在約化速度較大一端。

(3)數(shù)值模擬的結(jié)果與線性化分析和實驗結(jié)果一致，并且表明線性化分析得出的鎖定區(qū)間范圍的大小與質(zhì)量比仍成負(fù)相關(guān)關(guān)系的結(jié)果可推廣至低質(zhì)量比情況下。此外，數(shù)值模擬不僅捕獲了渦激振動中的“拍”現(xiàn)象，并且可反應(yīng)出上下端分支位移與升力之間相位發(fā)生180°變化。