袁方正， 袁德奎， 曾 攀， 劉長(zhǎng)根

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

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

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

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

1 線性化分析

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

圖1 單自由度橫向振動(dòng)的物理模型

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

(1)

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

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

y=Asin(ωext)。

(2)

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

CL=Clsin(ωext+φ)。

(3)

式中：Cl為升力系數(shù)的幅值；φ表示升力系數(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)。

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

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

(7)

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

(8)

(9)

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

2 實(shí)驗(yàn)結(jié)果分析

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

圖2 Williamson實(shí)驗(yàn)結(jié)果

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

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

3.1 流動(dòng)控制方程

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

(10)

3.2 湍流模型

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

(11)

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

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

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

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

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

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

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

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

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

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

(12)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.3 組合因子對(duì)渦激振動(dòng)特性的影響

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

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

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

5 結(jié)論

本文在前人實(shí)驗(yàn)結(jié)果的基礎(chǔ)上，采用線性化分析和數(shù)值模擬相結(jié)合的方法，對(duì)影響渦激振動(dòng)特性的參數(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)時(shí)，其振幅大小由組合因子m*ζ決定。組合因子越小，圓柱的振幅越大，反之則越小，但是組合因子對(duì)振幅的影響幅度有限；柱體鎖定區(qū)間范圍的大小由質(zhì)量比m*決定，質(zhì)量比越小，鎖定范圍越大，反之則越小，質(zhì)量比對(duì)鎖定區(qū)間范圍的大小有顯著影響，影響的主要范圍是在約化速度較大一端。

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