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

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

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

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

本文在已有研究的基礎上，對流體荷載進行線性分解推導得到了鎖定區(qū)間柱體響應的頻率比、振幅與流體荷載和系統(tǒng)物性參數(shù)之間的定性關系。根據(jù)推導得出的結(jié)果，結(jié)合相關的經(jīng)典實驗，對圓柱渦激振動特性的影響因素(質(zhì)量比，阻尼比，質(zhì)量阻尼比)開展了系統(tǒng)的數(shù)值實驗，對比了不同因素下圓柱鎖定區(qū)間、振幅、升阻力系數(shù)等關鍵參數(shù)和特征，進一步分析了不同質(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ū)，將圓柱的振動近似認為簡諧振動，且泄渦頻率完全被固定在圓柱振動頻率上。因此假設圓柱橫向位移y的表達式為：

y=Asin(ωext)。

(2)

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

CL=Clsin(ωext+φ)。

(3)

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

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ì)量有關。

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

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

(7)

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

(8)

(9)

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

2 實驗結(jié)果分析

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

圖2 Williamson實驗結(jié)果

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

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

3.1 流動控制方程

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

(10)

3.2 湍流模型

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

(11)

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

3.3 模型設置及求解方法

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

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

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

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

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

首先分析了邊界網(wǎng)格(主要是第一層的網(wǎng)格無量綱高度y+和圓周網(wǎng)格數(shù))對計算精度和效率的的影響，并從中選擇最接近實驗結(jié)果的壁面邊界處理方法。參數(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)格第一層高度進行估算，得到結(jié)果后再不斷調(diào)整網(wǎng)格直到達到符合條件的y+值。本文設置了兩組工況分別研究了y+、圓周網(wǎng)格數(shù)不同引起結(jié)果的差異，計算結(jié)果見表1。

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

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

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

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

為了對渦激振動的影響因素進行分析，首先研究了固定阻尼比ζ=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ì)量比圓柱的振幅響應

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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