兩自由度運動圓柱繞流的離散渦方法模擬

董 婧，宗 智，李章銳，孫 雷，陳 偉

（1大連理工大學 工業裝備結構分析國家重點實驗室 船舶與海洋工程學院，遼寧 大連116024；

2中國艦船研究設計中心，武漢430050）

1 引 言

將鈍體放置在具有一定速度的來流當中，流體會繞過結構物，在其后面形成尾流區,并在一定條件下從結構兩側交替地發生周期性泄渦現象。而旋渦的產生和脫落的作用，使物體受到橫向和流向的脈動壓力。對于彈性圓柱體或彈性支承的剛性圓柱體，所產生的脈動力會引起結構物的振動，反過來結構物的振動也會影響到流場，這種流固耦合現象稱為渦激振動（Vortex Induced Vibration）[1]。當渦泄頻率接近于結構的某一固有頻率時，旋渦的泄放過程被結構的振動所控制,而使旋渦泄放和結構振動具有相同的頻率，并在較大的流速范圍內保持不變，這種現象稱為“鎖定”[2-3]。在發生鎖定時物體的振幅顯著增加，相位角發生突變。鎖定現象會在旋渦脫落頻率的很大范圍內發生，鎖定現象的機理尚不明確，大量的實驗數據表明，通常當折合速度（Ur）在4.5到8之間時發生鎖定現象。

海洋工程中的水下結構物經常采用圓柱體，其在波流作用下經常會發生渦激振動。當泄渦頻率接近于結構的某一固有頻率時，將會出現垂直于水流和圓柱體軸線方向的強烈振動，并且導致沿水流方向的振動和曳力增大。結構在兩個方向上的強烈振動將引起震蕩應力[4]。如果這種振動處于定常狀態而連續發生，將會導致結構的疲勞破壞，降低其使用壽命。渦激振動是海洋立管遭受破壞的主要原因，另外，在水流沖刷及海底表面不平等環境因素作用下，鋪設在海底的油氣管道難免形成懸空段，直接暴露在水流作用下也會發生渦激振動。

關于結構渦激振動的研究，國外起于20世紀60年代，而我國則從80年代開始。研究范圍涉及風工程和海洋工程；研究方法在早期側重于試驗研究，經過試驗研究和理論分析并重的階段，目前有重視理論分析的傾向[2-4]。眾所周知，在通常的渦激振動中橫向振動占主導地位，振幅大大超過流向振動，所以在研究過程中，很多學者采用簡化的單自由度運動模型，只考慮橫向振動而忽略了流向振動因素。那么這種做法是否科學合理？流向振動到底對渦激振動有何影響？本文的目的之一就要驗證它。此外，本文還要對渦激振動的特征機理進行分析探討。

本文采用離散渦方法（Discrete Vortex Method，DVM）來研究彈性支承的二維圓柱繞流的渦激振動（VIV）問題。對于圓柱的運動情況的研究分為單向自由度的橫向運動和兩向自由度的流向和橫向共同運動兩種情況。本文分別進行了單向自由度橫向運動和兩向自由度橫向和流向共同運動這兩種情況的模擬，并討論了兩者之間的聯系，在一系列共16組（單向、兩向運動情況各8組）數值計算中采用同Stappenbelt[6]的試驗一致的各項參數，其中包括相同尺度的圓柱模型，均為單位長度，直徑0.055 4m；從2.36到12.96變化的8個不同的質量比，以及隨質量比的改變而改變的水中固有頻率；從3到16變化的折合速度，以及隨折合速度和固有頻率的改變而改變的速度和雷諾數；還有不變的阻尼比，按試驗值均為0.006。

本文首先對受力系數和響應位移進行了時域分析。從脈動升力和圓柱橫向響應的振幅、周期和相位關系來判斷鎖定區域和非鎖定區域，討論在不同折合速度下，振動幅值的初始分支、上端分支和下端分支的變化規律，以及質量比的增加對振動幅值、鎖定位置和鎖定區域的影響。

接下來進行的是頻域響應的分析，就泄渦頻率和振動頻率關系與Stappenbelt實驗結果進行了對比，吻合完好，斯托哈爾數都處于穩定范圍內，并且進一步討論了質量比的增加對振動幅值、鎖定位置和鎖定區域的影響。

最后進行了不同折合速度、不同自由度下尾渦模式的分析，并結合一些文獻中已發表的結果進行了對比和討論。在分析的過程中認為數值模擬的絕大部分尾渦模式與實驗和計算相一致，并且發現了新的尾渦模式。

2 離散渦方法的原理[1，5]

假設流體為不可壓縮流體，那么它的控制方程為N-S方程：

因渦量ω=▽×V，對上式取旋度，得：

如只研究二維平面流，則ω·▽V=0，由此得：

離散渦方法是將連續的渦量場離散，在拉格朗日框架下進行求解，著眼于追蹤每個渦元粒子運動的時間歷程。對于不可壓縮的粘性流動，上世紀90年代Chorin[7]提出了算子分裂法，渦量方程（2）可分裂成對流和擴散兩部分分別進行求解。

對流部分：

粘性擴散部分：

其中，υ為流體的運動粘性系數，ω為渦量，ψ為流函數。

方程（4）代表二維不可壓縮無粘性流動，如果當無窮遠處的自由流動條件已知，流場中不存在外部邊界時，它的解就是著名的Biot-Savart定律：

式中 r、r′，表示矢徑，V∞代表無窮遠處的速度。

而擴散部分與粒子的布朗運動的解具有相同的形式，在無粘解的基礎上，可以采用隨機走位的方法來模擬，下面具體介紹隨機走位的方法。

方程（5）為一維熱擴散型方程，它的基本解是格林函數：

這與一個均方差為零、標準差為σx的隨機變量的概率密度函數在形式上是一致的，即：

其中N是流場中渦元總數，下標i為渦元編號，Pi和Qi分別是在（0,1）和（0，2π）區間內均勻分布的兩個相互獨立的隨機數。

假設t時刻渦元的位置xi（t），yi（t）已給定，則t+Δt時刻渦元的隨機走位為

在此應當指出的是，每個渦元的環量Γi是不隨時間變化的。

3 離散渦方法的數值實現

離散渦方法的實質是用離散的點渦來模擬物體的邊界層和剪切層，認為粘性影響集中在物體附近很薄的邊界層以及分離后的剪切層內，其渦量來自于不斷分離出來的點渦供給，在整個流場中渦量守恒。離散渦方法主要分為物面渦的生成、渦元的對流和擴散、渦元的消除和融合，物體受力計算等幾個方面。

首先，將連續的渦量離散成為若干個渦元（如圖1所示），根據流場條件計算出每個渦元的位置和強度。

其次，要求出流場中的速度分布，每個渦元的速度都由兩部分疊加而成，一部分為定常流場的繞流速度；另一部分為渦元之間的誘導速度，可以通過泊松方程▽2ψ=-ω和邊界條件聯立求解。

在求解邊界條件時，本文采用流函數列式的方法即渦元法（Vortex Element Method）[1，5]。 渦元法利用物面是一條流線

圖1 離散渦元示意圖Fig.1 Position of vortex generation

來滿足物面條件，通過離散的流函數列式，可以確定各新生渦元的強度。Spalart等人在1983年引入了剛性渦，其半徑為σ，周向的誘導速度為：

當渦元的誘導速度為（5）式時，其相應的流函數為

則在控制點k處的流函數為：

其中N和M分別為新生渦元總數和尾渦渦元總數。

當圓柱靜止不動時物面為一條流線，則

將（13）式代入（14）式可以得到一組以環量Γ為變量的線性方程組：

其中系數矩陣AI僅與控制點和新生渦元的位置有關，列向量b由來流和尾流中的渦元分布來決定。

關于物體的受力，Quartepelle和Napolitanp在1983年推導了用整個流場的信息來表示的受力系數公式[9]：

對于運動的圓柱，（14）式應改寫為：

其中VB圓柱運動的速度矢量，n為控制點k處的物面法向。相應的（15）式中的列向量b中也會增加與速度有關的因素，而系數矩陣AI保持不變。此時（16）式必須加入圓柱的慣性力項，于是改寫為：

無論圓柱運動與否，現有渦元的速度都由環量Γ表示如下：從上式中可以看出，在求取每個渦元的速度時，都需要計算其他所有渦元對它的誘導速度，如果把計算單個渦元速度的運算次數計為NV，那么在一個時間步內計算所有渦元的速度，則需要N2V次。而每一時間步內都有固定數目的新生渦元進入尾流中，隨著時間的增長NV直線增加，其結果必然導致單個時間步長的計算時間越來越長，計算緩慢。因此，必須有效控制尾流中渦元的數目。同時，在渦的對流和擴散中，渦元會越過邊界移動到圓柱內部引起混亂，也必須予以消除。Spalart提出的渦融合機制認為，如果兩渦元滿足下面（20）式的條件就進行合并。

其中z=x+yi為渦元在復平面內的位置，D0和V0都是控制參數，其中D0∝U∞，控制圓柱附近的渦元數目；V0∝10-6U∞，控制流場中渦元總數。合并后新渦元的位置為（21）式，強度為（22）式。

運動圓柱的渦激振動問題可以從兩方面入手，即強迫振動和自激振動。本文進行了自激振動的模擬，分別采用單向自由度和兩向自由度的彈簧-阻尼質量系統來模擬二維圓柱，如圖2所示。

圓柱在橫向和流向的運動方程可分別表示為：

圖2 彈性支承的圓柱模型Fig.2 Vibration model of elastic mounted cylinder

其中m為圓柱的質量，c為系統阻尼系數，k為彈簧剛度系數。Fx（t），Fy（t）分別為流向和橫向的流體力，它們可以分別用阻力系數和升力系數表示為：

將（23）、（24）式轉化成無量綱形式為：

采用4階精度的Runge-Kutta法進行迭代求解。

4 數值結果和實驗對比

在兩向自由度橫向振動的數值模擬中，為方便與其它試驗和數值結果對比，本文選用Stappenbelt[3]的實驗參數，如表1所示

除表1所列的參數之外還有其他一些在這8組實驗當中不變的參數，包括圓柱直徑D為0.055 4m，折合速度Ur從3變化到16，以0.5遞增，時間步長Δt=0.1D/U，每時間步有新生渦120個，總共1 200步。雷諾數Re隨折合速度Ur變化，其值在2.84×105到1.52×106之間，每組實驗都分為單自由度橫向振動和兩自由度流向和橫向耦合振動兩種情況計算，為方便表述，分別記為1dof和2dof。

在圖 3（a）到（f）給出了質量比為 2.36時的 6個典型時刻的升力系數CL、阻力系數CD、無因次化的橫向位移Y/D、流向位移X/D的時程曲線。通觀這6幅圖可以得到以下結論：

表1 振動模型計算參數Tab.1 Parameter of vibration model

從受力系數來看，阻力系數CD始終處于一個較高的水平，其波動的幅值較小，這其中穩定的部分是由均勻來流的持續沖擊所致，而波動的部分是由泄渦產生的脈動壓力所致，所以從上面的幾幅圖的比較可以看出，折合速度越大，阻力系數的平均值就越高，這是符合常理的。升力系數CL的平衡位置在0附近，而波動的幅值相對較大，這是由于升力系數只與泄渦產生的脈動壓力有關，橫向的脈動壓力幅值大于流向的脈動壓力幅值，而均勻來流的持續沖擊力對升力系數沒有直接的影響。升力系數CL的脈動周期始終約為阻力系數CD脈動周期的2倍，這也是由泄渦產生的脈動壓力所致。從位移曲線來看，流向位移X/D受阻力系數CD影響，在一個較高的水平上波動，在第一個周期里有較高的幅值，這是由于在靜止狀態突然加力所致，在后面的幾個周期里就會平穩下來，不過這也在另一個側面說明了本文計算程序的穩定性。橫向位移Y/D受升力系數CL影響，則在0軸上下波動。

圖3（a）中曲線對應折合速度Ur=4.5，處于初始分支。此時升力系數CL和阻力系數CD出現階段性的脈動，周期大概在25左右，這與彈簧—阻尼系統在兩個方向的固有頻率接近有關。橫向位移Y/D與升力系數CL的周期和相位都是一致的，這一現象在初始分支中普遍存在，并與眾多試驗和數值模擬結論相一致[10-12]。圖3（b）中曲線對應折合速度Ur=6，處于初始分支和上端分支的臨界狀態，橫向位移Y/D的振動幅值較大且振動平穩，但相位和周期仍然符合初始分支的特點。圖3（c）中曲線對應折合速度Ur=7.5，處于上端分支，橫向位移Y/D的振動幅值較大且振動平穩，橫向位移Y/D與升力系數CL的周期一致而相位角相差180°，這是上端分支鎖定區域的典型特點。圖3（d）中曲線對應折合速度Ur=9，處于上端分支和下端分支的臨界狀態，此時橫向位移Y/D與升力系數CL仍然是周期一致而相位角相差180°，而橫向位移Y/D的振動幅值開始下降且振動不穩定。圖3（e）、（f）處于下端分支的中段和末段，在這一階段，橫向振動逐漸減退，流向振動和橫向振動的幅值相當。

圖3 升力系數CL、阻力系數CD和流向位移X/D、橫向位移Y/D的時程曲線，其中質量比m*=2.36，水中固有頻率Fn=1.711Fig.3 History of lift coefficient CL,drag coefficient CD,stream-wise displacement X/D,and transverse displacement Y/D,here mass-ratio m*=2.36,natural frequency in water Fn=1.711

圖4到圖7為不同質量比時單自由度1dof和兩自由度2dof的橫向運動無因次振幅A/D同實驗值的對比，橫坐標為折合速度Ur。為與實驗的描述一致，這里用T表示橫向運動，T&S表示橫向和流向耦合運動。從圖中可以看出，本文計算結果與Stappenbelt實驗值的總體趨勢向一致，會出現明顯的初始分支、上端分支和下端分支。在兩自由度的振動中質量比較低時還會出現超上端分支，而后振幅驟降。從相同質量比時1dof和2dof振動情況的比較來看，相同折合速度下2dof振動時的幅值更大，這一特征隨著質量比的減小而增加，這說明了流向振動的參與對橫向振動會有促進作用。而從不同質量比之間的對比可以看出，隨著質量比的增加，無論是1dof還是2dof振動，其幅值都會隨之下降。本文計算結果的不足之處在于模擬到的下端分支處的振幅偏高。

圖4 m*=2.36時單雙自由度橫向無因次振幅A/D同實驗對比Fig.4 Non-dimensional transverse amplitudes A/D of 1dof and 2dof compared with experimental results for m*=2.36

圖5 m*=3.68時單雙自由度橫向無因次振幅A/D同實驗對比Fig.5 Non-dimensional transverse amplitudes A/D of 1dof and 2dof compared with experimental results for m*=3.68

圖6 m*=5.19時單雙自由度橫向無因次振幅A/D同實驗對比Fig.6 Non-dimensional transverse amplitudes A/D of 1dof and 2dof compared with experimental results for m*=5.19

圖7 m*=6.54時單雙自由度橫向無因次振幅A/D同實驗對比Fig.7 Non-dimensional transverse amplitudes A/D of 1dof and 2dof compared with experimental results for m*=6.54

圖8共有10組質量比的單自由度運動時無因次振動頻率和泄渦頻率圖，橫坐標為折合速度。其中前兩幅圖為Stappenbelt實驗中已發表的數據。后8幅圖為本文計算結果，與實驗結果的趨勢一致，在鎖定區域以內，泄渦頻率fs和圓柱振動頻率fv一起鎖定在固有頻率fn附近；在鎖定區域以外，泄渦頻率fs沿折合速度會形成一條直線，在這條直線上有一個相對穩定的斯托哈爾數St，這將在表2中說明。

表2為單自由度振動模型的鎖定參數 （包括初始鎖定的折合速度Ur和鎖定區域折合速度的范圍ΔUr）及斯托哈爾數St的計算結果同Stappenbelt實驗值的對比。從折合速度Ur的比較中可以看出，兩組數據的大小相接近，且都隨著質量比的增加而增加。從鎖定區域的范圍ΔUr可以看出，本文計算與實驗結果符合得較好，最多相差0.6，且都隨著質量比的增加而減小。從斯托哈爾數St的比較可以看出，本文計算與實驗結果符合得較好，最多相差0.08，且都隨著質量比的增加而增大。雙自由度振動模型的計算結果與之類似。

圖8 無因次振動頻率和泄渦頻率關系圖Fig.8 Non-dimensional vibration and vortex shedding frequencies

表2 鎖定及斯托哈爾數計算結果同實驗對比Tab.2 Lock-in values and Strouhal number compared with experimental results

橫向振動響應分支在振幅上有較大變化，相應的尾渦模式也有明顯的差異。在本文的計算結果當中，Re在105以上，無論是1dof還是2dof的振動，在初始分支中都表現出2S模式（two single vortex shedding per cycle）如圖9所示，而通常情況下在上端分支和下端分支所對應的尾渦則為2P模式（Two pairs vortex shedding per cycle）。上端分支所對應的尾渦中2P模式為一強一弱的兩個旋渦如圖10所示，其中小的旋渦會在尾流中不斷的對流和擴散中消失，在流場中較遠處只剩下較大的旋渦；而下端分支所對應的尾渦中2P模式為強度相當的兩個旋渦，如圖11所示，這兩個漩渦在流場較大范圍內運動都不會消失。這與黃志勇[11]文獻中描述的結果類似。而本文的計算中還發現了，在個別上端分支中在一個周期內先是產生兩條細長的渦帶而后分化成4個小的旋渦，暫時把這種尾渦模式也歸為2P模式的一種，如圖12所示，本文認為這是從2P模式到2T模式的一種過渡模式，其生成機理有待于進一步研究。

圖9 初始分支的2S瀉渦模式Fig.9 2S model of vortex shedding in the initial branch

圖10 上端分支中兩旋渦不等的2P瀉渦模式Fig.10 2P model of vortex shedding in the upper branch(each pair of vortex blobs with unequal scales)

圖11 下端分支中兩旋渦相等的2P瀉渦模式Fig.11 2P model of vortex shedding in the lower branch(each pair of vortex blobs with equal scales)

圖12 特殊的2P瀉渦模式Fig.12 The special 2P model of vortex shedding

圖13 超上端分支中的2T瀉渦模式Fig.13 2T model of vortex shedding in the super upper branch

2T模式（Two triplets of vortices per cycle）在一個泄渦周期內產生兩組渦團，每個渦團產生三個旋渦，其中兩個旋渦的旋向與另一個相反。出現在兩自由度振動的上端分支的頂點處，也就是所說的超上端分支，如圖13所示。

在一些過渡階段還會出現以上多種尾渦模態交替出現的情況，是由于運動圓柱的位移和速度的變化所致，在這里不做詳細介紹。

5 結 論

本文利用離散渦方法對二維彈性支承圓柱繞流渦激振動問題進行了較為詳細的計算，可以得到以下結論：

首先對受力系數和響應位移進行了時域分析。從脈動升力和圓柱橫向響應的周期和相位關系可以看出有明顯的鎖定區域和非鎖定區域，觀察到在不同折合速度下，振動幅值有著明顯的初始分支、上端分支和下端分支的變化，并且隨著質量比的增加，振動幅值會有所下降，這與眾多觀察和實驗的本質特征一致。

在單自由度和兩自由度振動幅值的對比中發現，兩自由度振動在上端分支處的幅值明顯高于單自由度的情況，說明流向振動對橫向振動有促進作用，用單自由度建立的渦激振動模型預報結果將會偏于保守，在此建議在建立渦激振動模型時，流向振動因素應當予以重視。

從進一步的頻域分析結果來看，泄渦頻率和振動頻率關系與Stappenbelt實驗結果吻合完好，斯托哈爾數都處于穩定范圍內，并且驗證了隨著質量比的增加，初始鎖定時的折合速度會提高，而鎖定區域的范圍會減小。

在進行尾渦模式分析時與黃志勇在文獻中已發表的結論基本一致，并且在上端分支中發現了新的尾渦模式—有4個小的旋渦的2P模式。本文認為這是從2P模式到2T模式的一種過渡模式。

[1]Lewis R I.Vortex element method for fluid dynamic analysis of engineering systems[M].Cambridge:Cambridge University Press,1991.

[2]Gabbai R D,Benaroya H.An overview of modeling and experiments of vortex-induced vibration of circular cylinders[J].Journal of Sound and Vibration,2005,282(3):575-616.

[3]Liao Jungchi.Vortex-induced vibration of slender structures in unsteady flow[M].Cambridge,USA:Massachusetts Institute of Technology,2002.

[4]Pan Zhiyuan,Cui Weicheng,Liu Yingzhong.A predicting model for self-excited VIV of a circular cylinder at low massdamping[J].Journal of Ship Mechanics,2005,10(5):115-124.(in Chinese)

[5]Saltara F.Meneghini J R,Fregonesi R A.Numerical simulation of flow around elastically mounted cylinder[J].International Journal of Offshore and Polar Engineering,2003,13(2):99-104.

[6]Stappenbelt B,Lalji F.Low mass ratio vortex-induced motion[C].16th Australasian Fluid Mechanics Conference,2007,12:1491-1497.

[7]Chorin A J.Numerical study of slightly viscous flow[J].Journal of Fluid Mechanics,1973,57:785-796.

[8]Lewis R I.Vortex element method for fluid dynamic analysis of engineering systems[M].Cambridge:Cambridge University Press,1991.

[9]Chorin A J.Numerical study of slightly viscous flow[J].Journal of Fluid Mechanics,1973,57:785-796.

[10]Huang Zhiyong,Pan Zhiyuan,Cui Weicheng.Numerical simulation of VIV of a circular cylinder with two degrees of freedom and low mass-ratio[J].Journal of Ship Mechanics,2007,11(1):1-9.(in Chinese)

[11]Khalak A,Williamson C H K.Dynamics of a hydroelastic cylinder with very low mass and damping[J].Journal of Fluids and Structures,1999,13:813-851.

[12]Zhou C Y,So R M C,Lam K.Vortex-induced vibrations of an elastic circular cylinder[J].Journal of Fluids and Structures,1999,13:165-189.