低雷諾數下直圓柱和波型圓柱受迫振動的數值研究

平 煥， 張 凱,2， 周 岱,3,4， 包 艷， 朱宏博， 韓兆龍

(1.上海交通大學 船舶海洋與建筑工程學院，上海 200240； 2.橫濱國立大學 都市創新學院，橫濱 2408501；3.高新船舶與深海開發裝備協同創新中心，上海 200240； 4.上海交通大學 海洋工程國家重點實驗室，上海 200240)

自從Bishop等[1]開創性地開展直圓柱的受迫振動實驗以來，均勻來流下直圓柱的受迫振動問題成了渦激振動領域一個聚焦的熱點[2-4]。相對于自激振動，受迫振動能更多地反映出尾渦與柱體運動之間的相互作用以及能量轉移的關系，因此具有重要的研究意義。

在受迫振動的研究中，發現了許多重要的現象，例如：Ongoren等[5]通過實驗研究發現，當直圓柱的振動頻率(fex)與自然瀉渦頻率(fs)之比fex/fs從0.85增至1.17時，初始瀉渦形成的相位會發生近180°的轉變，即在圓柱一側形成的渦在圓柱運動到了同一側的極限位置處脫落的現象轉變為在圓柱運動到了另一側極限位置處才脫落。Gu等[6]研究了直圓柱在Re=185和Re=5 000下，振幅為0.2D(D為直圓柱直徑)的受迫振動，也發現了類似的現象；受迫振動下，直圓柱尾流區并不是單一的瀉渦模式，Williamson等[7]通過大量研究統計，觀察到了2S，2P和P+S三種主要的瀉渦模式，并在A*-fex/fs平面上進行了劃分，同時給出了臨界曲線。

另外一個備受關注的現象就是“鎖定(Lock-in)”，即當運動圓柱的振動頻率接近自然瀉渦頻率時，圓柱后側的瀉渦頻率(fcl)不再由斯特勞哈爾數St確定，而是鎖定到振動頻率上。但對于瀉渦頻率鎖定到振動頻率的具體判定，不同的學者有不同的定義，例如：Karniadakis等[8]研究了Re=100下直圓柱的受迫振動。他們將尾流區各個監測點上的流向速度能量譜的主頻等于振動頻率的情況視之為鎖定；Meneghini等則將升力系數能量譜的主頻等于振動頻率定義為鎖定；Kumar等[9]提出鎖定必須同時滿足：①升力系數能量譜的主頻等于fex；②主頻以外的其他頻率等于fex的整數倍這兩個條件；另外，Koopmann[10]發現振動頻率在一定范圍時，直圓柱后側由傾斜瀉渦變為平行排列瀉渦時，就說明發生了鎖定。不同的定義對數據結果的呈現產生了一定的差異，特別是對鎖定區間的劃分。

近些年研究發現，波型圓柱(直圓柱橫截面直徑沿軸向按正弦曲線形式變化得到的一種圓柱體)具有較好的減阻和抑制升力的效果，展現出了潛在的工程應用價值，故開始引起學者們的關注。Ahmed等[11]通過實驗研究了不同波長波型圓柱的表面壓力分布，發現最大截面處的壓力系數要大于最小截面處。鄒琳等[12]通過數值模擬研究了Re=3 000下不同表面波動幅值的波型圓柱的物理特性，指出較小的表面波幅并不能達到減阻的目的。但隨波幅的增加，平均阻力系數開始減小，最大減阻量可達16%。Lam等[13]則系統地研究了Re=100時不同波長和波幅的靜止波型圓柱繞流，發現在最優波長和波幅下卡門渦街能夠完全得到抑制，同時定義了A，B，C三種瀉渦模式。

綜合以上國內外研究現狀可以看出，目前受迫振動的研究對象都是直圓柱，幾乎沒有關于波型圓柱受迫振動的研究，同時對兩者在受迫振動下物理特性比較的文獻也較稀少。因此，本文在低雷諾數下，對均勻來流中作橫向正弦受迫振動的直圓柱和波型圓柱進行了數值模擬。以Re=150為例，對比研究直圓柱和波型圓柱在不同振幅和振動頻率下的阻、升力和尾渦的變化特征以及鎖定區間的確定。

1 數值模型

1.1 幾何模型

圖1給出了波型圓柱的模型示意圖，其直徑沿軸向Z的變化由下式給出

Dz=Dm+2acos(2πz/λ)

(1)

Dm=(Dmax+Dmin)/2

(2)

式中：Dz為波型圓柱沿軸向的局部直徑，Dm為平均直徑，Dmax，Dmin分為沿軸向最大和最小截面所在處的直徑，a為波型圓柱表面正弦曲線的幅值，λ為波長。由于a=0時，Dz=Dm=D，波型圓柱即成為直圓柱，故這部分只以波型圓柱為例作為說明。本文采用λ/Dm=2，a/Dm=0.18的波型圓柱研究其在受迫振動下的特性。

圖1 波型圓柱的模型示意圖

1.2 流動控制方程

三維不可壓縮均勻黏性流體運動的基本控制方程為連續性方程和Navier-Stokes方程，在直角坐標系下，分別可以表示為如下的無量綱形式

(3)

(4)

式中：下標i為坐標分量(三維情況下，i=1，2，3為x，y和z三個坐標分量)，ui為i方向的流體運動速度分量，p為壓力，Re=U∞Dm/υ為雷諾數，U∞為均勻來流速度，υ為流體的運動學黏性系數。對于圓柱運動引起的網格變形問題則通過求解Laplace方程來解決

Δ·(γ▽dm)=0

(5)

式中:dm為網格位移，γ為到圓柱表面的逆距離控制的擴散系數。本文利用CFD軟件OpenFOAM內置求解器pimpleDyMFoam求解上述方程。具體地說，采用有限體積法求解N-S方程：對時間項采用隱式積分方法；對流項采用二階迎風離散格式；對控制方程中的速度和壓力耦合則是采用PIMPLE算法來處理。PIMPLE算法是SIMPLE和PISO算法的結合體，有關該算法的具體介紹，請參閱相關文獻[14-15]。

1.3 計算區域及邊界條件

1.3.1 計算區域

采用圓柱體型計算區域，坐標原點位于圓柱底面的中心，計算區域半徑為30Dm，高度為2Dm。Lam等[16]的研究表明，無論流動是處于層流還是湍流，選取一個波長的計算區域高度可以得到準確的模擬結果。另外，本文對最大振動頻率下的工況(A*,f*)=(0.25，0.3)進行了驗證，發現采用兩個波長的計算區域高度得到的結果與采用一個波長時一致?；谏鲜隼碛?，本文選取2Dm(一個波長)的計算區域高度能夠保證獲得準確的計算結果。沿x-y平面，網格呈非均勻分布，靠近圓柱表面網格密集，遠離圓柱，網格逐漸變稀疏。而沿圓柱的軸向Z，網格均勻分布。波型圓柱周圍的網格分布，如圖2所示。

圖2 波型圓柱周圍的網格分布

1.3.2 邊界條件

1.4 無量綱物理量的定義

本文中用到的一些重要的無量綱物理量定義，如表1所示。

2 網格及算法可靠性驗證

對直圓柱分別采用較為稀疏的網格模型mesh I、中等疏密的網格模型mesh II和較為稠密的網格模型mesh III進行網格可靠性驗證，相應的網格參數，見表2。

表3給出了在Re=100下靜止直圓柱繞流計算得到的升力系數均方根值CL(rms)、平均阻力系數CD(mean)、斯特勞哈爾數St以及與其他已有參考文獻結果的比較?？梢园l現，使用網格模型mesh II和III得到的計算結果與文獻結果較為接近，說明這兩種網格模型均滿足計算精度的要求?？紤]到計算效率，對于直圓柱，選用網格模型mesh II作為后續計算的網格模型。

表1 無量綱物理量的定義

表2 直圓柱網格模型參數

表3 直圓柱網格模型計算得到的結果和已有文獻結果的比較

為進一步驗證算法的正確性，本文對直圓柱作橫向受迫振動的情況進行了驗證。研究了Re=185，A*=0.2下6種不同振動頻率的工況。圖3為本文數值結果與Guilmineau等[20]結果的對比，其中CD(rms)為阻力系數均方根值?？梢钥闯?，兩者結果相吻合且具有一致的變化趨勢，說明本文采用的算法適用于低雷諾數下的受迫振動。

另外，我們對波型圓柱計算網格模型的可靠性進行了驗證。同樣采用三種不同疏密程度的網格模型mesh IV、mesh V、mesh VI，相應的網格參數見表4。

圖4則給出了波型圓柱在(A*,f*)=(0.25, 0.15)時，通過三種網格計算得到的阻、升力系數的時歷曲線?？梢钥吹?，三條時歷曲線基本一致。計算表明，由mesh IV獲得的CL(rms)、CD(rms)和CD(mean)相對于mesh VI，其誤差分別為4.9%、1.9%和0.6%；而由mesh V獲得的結果相對于mesh VI，誤差則分別為1.9%、1.8%和0.2%，說明計算結果在三種網格上趨于網格收斂?？紤]到計算效率，對于波型圓柱，選用網格模型mesh V作為后續計算的網格模型。

圖3 受迫振動下直圓柱計算得到的結果和已有文獻結果的比較

Fig.3 Comparison with published results for uniform circular cylinder undergoing forced oscillation

表4 波型圓柱網格模型參數

圖4 不同波型圓柱網格模型計算得到的阻、升力時歷曲線

Fig.4 Time histories of drag and lift coefficients with three different grid resolutions for wavy cylinder

3 計算結果及分析

3.1 靜止波型圓柱繞流的數值模擬

圖5為Re=150時，模擬靜止波型圓柱繞流得到的阻、升力系數時歷曲線。從圖5可知，阻力系數無波動，升力系數為零，說明渦脫落現象被完全抑制，卡門渦街型瀉渦頻率消失。圖6則為計算穩定后某一時刻的三維渦結構圖。從圖6可知，兩側對稱分布的自由剪切層穩定地朝下游延伸發展，沒有卷起形成渦對。故本文采用的波型圓柱在靜止狀態下能完全抑制瀉渦的發生。

圖5 靜止波型圓柱繞流的阻、升力時歷曲線

圖6 靜止波型圓柱繞流的三維渦結構圖(ωz)

Fig.6 Instantaneous vortical structure for flow past a stationary wavy cylinder

3.2 受迫振動下直圓柱和波型圓柱繞流的數值模擬

3.2.1 阻、升力系數

圖7給出了不同振幅下，直圓柱和波型圓柱的CL(rms)，CD(rms)以及CD(mean)隨振動頻率的變化情況。對于直圓柱，研究了區間為 [0.05,0.5]的10種不同振幅；對于波型圓柱，研究了A*=0.05, 0.25, 0.5三種有代表性的振幅，分別對應著小振幅、中振幅和大振幅的情況。

從圖7(a)可知，在f*<0.18的低頻范圍內，隨著振動頻率的增加，直圓柱的升力系數均方根值曲線先下降后上升，但不同振幅的曲線開始下降的點所對應的振動頻率不同，其隨振幅的增加而減小。通過對比圖8可知，該振動頻率即為該振幅下系統開始進入鎖定狀態的頻率，即升力系數均方根值的下降是由系統開始進入鎖定狀態引起的。在f*>0.18的高頻范圍內，升力系數均方根值對振幅的依賴性較大，振幅越大引起的升力系數幅值也越大。另外，值得注意的是，在振動頻率區間[0.18,0.2]內，各個振幅的曲線有一個突變，通過對比圖8可知，這是由于系統脫離鎖定狀態引起的。對于波型圓柱，從圖7(d)可知，A*=0.5與A*=0.25時的曲線走勢相似。初始階段由于振動頻率小，振動十分緩慢，接近于靜止狀態，而靜止狀態下波型圓柱的瀉渦是被完全抑制的，故升力系數均方根值不變且較小，升力是由圓柱的受迫振動引起。接著，由于瀉渦的發生，升力系數均方根值出現了小的陡增。隨著振動頻率的進一步增加，系統開始進入鎖定狀態(對比圖9)，導致曲線有一個明顯的下降。最后，伴隨著系統脫離鎖定狀態，升力系數均方根值出現突增。A*=0.05時的情況也類似，只是曲線對應階段的變化幅度比前兩者小。但有所不同的是，在f*=0.25, 0.26和0.3處波型圓柱后側的瀉渦被完全抑制，導致其升力值明顯小于相鄰值。對比CL(rms)縱坐標值可以發現，波型圓柱較直圓柱在低頻段體現出了抑制升力的效果。

振幅和振動頻率對阻力系數均方根值的影響示于圖7(b)和圖7(e)。對于直圓柱，總體上，振幅越大，平均阻力系數也越大。每條曲線經歷了兩次下降過程。第一次下降代表系統進入鎖定狀態；第二次下降代表系統脫離鎖定狀態；從圖7(e)可知，波型圓柱的情況與直圓柱類似。A*=0.5與A*=0.25時，曲線的兩次下降過程同樣分別對應著系統進入和脫離鎖定狀態。A*=0.05時的曲線在高頻段的表現則略有差異：隨著振動頻率的增大，曲線并沒有呈現上升趨勢，反而由于在f*=0.25，0.26和0.3處由于瀉渦被抑制，導致出現下凹。對比CD(rms)縱坐標值可以發現，波型圓柱在高頻段體現出了減阻的效果。

圖7(c)和圖7(f)則展示了平均阻力系數與振幅和振動頻率之間的關系。圖7(c)曲線的總體走勢和圖7(b)大致相似，但初始段不同振幅對平均阻力系數幾乎沒有影響；對于波型圓柱，通過對比圖7(f)和圖7(e)可知：平均阻力系數曲線少了第一次下降過程。

(a) 直圓柱CL(rms)

(b) 直圓柱CD(rms)

(c) 直圓柱CD(mean)

(d) 波型圓柱CL(rms)

(e) 波型圓柱CD(rms)

(f) 波型圓柱CD(mean)

3.2.2 鎖定現象

采用Kumar等對鎖定的定義，即：①升力系數能量譜的主頻等于fex；②主頻以外的其他頻率等于fex的整數倍。不同時滿足上述兩個條件的情況則視為非鎖定。

圖8給出了直圓柱的鎖定區間和A*-f*平面上4個特定工況的升力系數時歷曲線及對應的能量譜。分界線下方區域代表非鎖定，上方區域代表鎖定。對于一個固定振動頻率，當振幅超過某個臨界值時，系統發生鎖定，臨界振幅隨振動頻率遠離直圓柱的自然瀉渦頻率而增加，這與Koopmann、Anagnostopoulos[21]的研究結果相似。

在f*<0.18的區間，分界線變化較平緩。圖8(a)給出了(A*,f*)=(0.05, 0.14)的升力系數時歷曲線及其對應的能量譜。從圖8(a)可知，升力并不是按規律的正弦形式波動，而是存在一個“差拍”現象。對應的能量譜顯示，主頻并不等于外加振動頻率，而是出現在自然瀉渦頻率附近，故處于非鎖定狀態。當振幅增加到0.25時，此時主頻已經固定到了外加振動頻率處，且其他頻率等于外加振動頻率的整數倍，這從圖8(b)可知，系統由之前的非鎖定狀態轉變到了鎖定狀態，這時升力系數時歷曲線也呈現出了規律的正弦波動形式。

在f*>0.18的區間，分界線較陡。圖8(c)和圖8(d)對應著(A*,f*)=(0.05, 0.19)和(A*,f*)=(0.25, 0.19)的工況，分別處于非鎖定和鎖定狀態。值得注意的是，由于(A*,f*)=(0.05, 0.19)的工況靠近過渡帶，故能量譜出現了兩個幅值相當接近的頻率(見圖8(c)局部放大圖)，但由于主頻與外加振動頻率并不重合，故仍處于非鎖定狀態。

(a) (A*, f*)=(0.05，0.14)

(b) (A*, f*)=(0.25，0.14)

(c) (A*, f*)=(0.05，0.19)

(d) (A*, f*)=(0.25，0.19)

圖9給出了波型圓柱在振幅為0.05、0.25和0.5時的鎖定區間以及A*-f*平面上4個特定工況的升力系數時歷曲線和對應的能量譜，同時對比了相同條件下直圓柱的結果。從圖9可知，在圖示頻率區間內, 大部分波型圓柱的鎖定區間被包絡在直圓柱的鎖定區間內。波型圓柱左側分界線的走勢和直圓柱相似，然而右側卻不同，在A*=0.25處有一個明顯的轉折。

(a) (A*, f*)=(0.05，0.15)

(b) (A*, f*)=(0.25，0.15)

(c) (A*, f*)=(0.05，0.18)

(d) (A*, f*)=(0.25，0.18)

圖9(a)、圖9(b)、圖9(c)、圖9(d)分別對應(A*,f*)=(0.05,0.15)，(0.25,0.15)，(0.05,0.18)和(0.25,0.18)四個特定工況。升力系數能量譜顯示，圖9(b)、圖9(d)處于鎖定狀態，圖9(a)、圖9(c)處于非鎖定狀態。與直圓柱的情況相同，在鎖定狀態下，波型圓柱升力系數時歷曲線為規律的正弦波動形式，而非鎖定狀態下則捕捉到了明顯的“差拍”現象。除了(A*,f*)=(0.25,0.15)的情況，波型圓柱較直圓柱起到了良好的抑制升力的效果，在(A*,f*)=(0.05,0.15)，(0.05,0.18)和(0.25,0.18)時升力系數均方根值分別降低了43.75%、49.82%和32.20%。

圖10、圖11分別為A*=0.25時，直圓柱和波型圓柱在各種振動頻率下的升力譜。其中縱坐標mag代表對升力系數的時程信號做快速傅里葉變換(FFT)后，獲得的以某種頻率成分振動的幅度特性。

圖10 A*=0.25時直圓柱在各種振動頻率下的升力譜

圖11 A*=0.25時波型圓柱在各種振動頻率下的升力譜

對于直圓柱，大部分范圍內，升力譜由外加振動頻率和卡門渦街型瀉渦頻率共同控制。但在振動頻率區間[0.14,0.19]內，代表卡門渦街型瀉渦頻率的豎線中斷消失，由振動頻率唯一控制，即發生鎖定；對于波型圓柱，同樣可以清楚地發現卡門渦街型瀉渦頻率的存在，這說明在靜止波型圓柱繞流中消失的卡門渦街型瀉渦頻率在受迫振動中又得以顯現。但不同于直圓柱，其發生了兩處中斷，分別在[0.14,0.18]和[0.05,0.08]區間上。前者代表著鎖定，而后者則是由于不瀉渦導致的。

3.2.3 尾渦特性

重點觀察了鎖定區間內圓柱尾流區的瀉渦模式及尾渦形態。通過統計發現，直圓柱后側的瀉渦模式由振動頻率控制，振幅則幾乎沒有影響。振動頻率較小時(低頻段)，瀉渦呈現2S模式；振動頻率較大時(高頻段)，則呈現C(2S)模式，這種瀉渦模式在Singh等[22]的研究中也有發現。以振幅A*=0.25為例，圖12給出了直圓柱處于平衡位置并朝y軸正方向運動時刻尾流區的瀉渦形式。

f*=0.14下，直圓柱尾流區呈現出典型的2S模式，即直圓柱在運動到最高點和最低點時各瀉放一個旋轉方向相反的渦，如圖12(a)所示。

當振動頻率增加到f*=0.19時,如圖12(b)所示，在近尾流區，同f*=0.14一樣，其呈現出典型的2S模式；但在離開直圓柱一定距離后，旋轉方向相同的渦發生合并，呈帶狀向下流延伸，即C(2S)模式。同時尾跡區的側向寬度隨離直圓柱距離的增加而增加。

(a) 2S瀉渦模式

(b) C(2S)瀉渦模式

Fig.12 Instantaneous vorticity contours showing two kinds of modes for uniform circular cylinder in lock-in

對于所有處于鎖定狀態的波型圓柱，其后側的尾渦呈現出了大致相同的特征：總體看來，和直圓柱瀉渦形式類似，但在軸向出現了輕微的扭曲，開始顯現出三維效應。最大和最小截面的瀉渦形態并沒有明顯的不同，均呈現出2S模式。以(A*,f*)=(0.25,0.18)為例，圖13分別給出了三維尾渦結構圖以及最大和最小截面處的瀉渦模式。這種尾渦呈現出來的特征和Lam等[13]提出的I類波徑比下的瀉渦模式A十分相似。

4 結 論

本文利用CFD軟件OpenFOAM對比研究了Re=150下直圓柱和波型圓柱的橫向受迫振動，特別是對波型圓柱受迫振動的研究具有創新意義。通過改變振幅和振動頻率，觀察阻、升力的變化特征，確定鎖定區間的范圍以及鎖定狀態下的尾渦特性，得到的主要結論如下：

(1) 總體上，不同振幅下的波型圓柱和直圓柱的升力系數均方根值、阻力系數均方根值以及平均阻力系數隨振動頻率變化的規律大致相同。波型圓柱較直圓柱體現出了一定的減阻和抑制升力的效果。減阻主要體現在高頻段，抑制升力則主要體現在低頻段。

(2)大部分波型圓柱的鎖定區間被包絡在直圓柱的鎖定區間內。在鎖定狀態下，兩者的升力系數時歷曲線均呈現出規律的正弦波動形式，而非鎖定狀態下，均捕捉到了明顯的“差拍”現象。

(a) 三維渦結構圖(ωz)

(b) 最大截面處的瀉渦模式

(c) 最小截面處的瀉渦模式

(3)在靜止波型圓柱繞流中消失的卡門渦街型瀉渦頻率又顯現在受迫振動中。

(4)鎖定區間內，在直圓柱后側觀察到了兩種瀉渦模式。瀉渦模式由振動頻率控制，與振幅幾乎無關。振動頻率較小時，呈現2S模式；振動頻率較大時，呈現C(2S)模式。而對于波型圓柱，在鎖定區間內只觀察到一種瀉渦模式，該模式和Lam等提出的I類波徑比下的瀉渦模式A十分相似。

(5)本文的工作為今后進一步開展均勻來流中波型圓柱渦激振動的研究作了鋪墊與參考。