振動半圓柱尾流中的二維擺動水翼推進性能研究

張 曦，蘇玉民，王兆立

（哈爾濱工程大學 水下智能機器人國防科技重點實驗室，哈爾濱 150001）

1 引 言

經過數千億年的進化，魚類等海洋生物具有在水中高效快速游動的能力。特別是在復雜流場中，魚類通過改變運動姿態來保持高效率的游動。已經有學者發現，魚類能夠利用魚體附近的旋渦，并從旋渦中吸收能量[1]。反觀諸如螺旋槳等傳統推進裝置，它們都不能很好地適應復雜流場。近些年來，已經有許多軀體-尾鰭模式（Body and Caudal Fin,BCF模式[2]）的仿生推進裝置用于小型水下航行器。美國麻省理工大學已經研制出利用旋渦提高推進性能的水下航行器VCUUV[3]。在VCUUV的基礎上，美國海軍研制出了世界上最快、最靈活的"金槍魚"潛水器，該潛水器在實驗過程中的優異表現已經引起了各國的重視。因此對BCF模式仿生推進裝置在復雜流場中，特別是在受到旋渦影響時的推進性能進行研究就顯得非常必要。

許多學者采用數值或實驗手段研究了旋渦對BCF模式仿生推進裝置的影響。在這些研究中，BCF模式仿生推進裝置通常被簡化為同時作升沉運動和搖擺運動的擺動水翼。Streitlien[4]在流場中布置點渦，應用勢流理論研究了點渦對二維擺動水翼水動力性能的影響。研究發現二維擺動水翼能夠從點渦中吸收能量從而提高推進效率。Gopalkrishnan[5]將大展弦比的擺動水翼置于卡門渦街中，從而找到了幾種旋渦和水翼渦相互作用的模式。Beal[6]也作了類似于Gopalkrishnan的實驗，實驗結果表明擺動水翼能夠吸收旋渦中的能量從而提高推力。Zhu等人[7]采用三維非線性面元法研究了擺動尾鰭在魚體脫落渦影響下的水動力性能，指出尾鰭可以利用魚體脫落渦提高推進性能。

本文二次開發了基于有限體積方法求解雷諾平均納維—斯托克斯（RANS）方程的CFD軟件FLUENT，計算了振動半圓柱后方二維擺動水翼的水動力性能。在求解過程中采用了動網格技術。分析了旋渦相互作用模式，探討了半圓柱渦對擺動水翼推進性能的影響。

2 數值計算模型

如圖1所示，半圓柱通過沿y軸升沉運動產生旋渦。半圓柱振動滿足的運動規律可以表示為：

式中：Ad為半圓柱升沉運動的振幅，f為半圓柱的運動頻率。

如圖1所示，二維擺動水翼既沿著y軸升沉又繞著水翼首緣點搖擺。因此水翼的擺動運動規律[8]可以表示為：

式中：Af為水翼升沉運動的振幅，θ0為水翼搖擺運動的振幅，f為水翼的運動頻率。為了便于計算和分析，取半圓柱的運動頻率和水翼的運動頻率相等。Φ0為水翼升沉運動和搖擺運動之間的相位差，在本文中 Φ0=-π/2。

圖1 半圓柱和水翼的周期性運動Fig.1 The cycle motion of an oscillating semicircular and a flapping hydrofoil

擺動水翼的推力系數Cx，側向力系數Cy和以水翼首緣點為作用點的力矩系數Cm由相應的推力Cx、側向力Cy和力矩Cm沿水翼表面積分[9]確定：

式中：C0表示水翼的特征弦長，ρ是水的密度，V0為無窮遠處來流速度。

擺動水翼的推進效率[10]表示為：

式中：Cxm為推力系數Cx在一個周期內的平均值。

輸入功率系數表示為：

3 數值計算方法

3.1 控制方程

對于不可壓縮的粘性流動，在笛卡爾坐標系下，用張量的形式表示的時均連續性方程和RANS方程可以寫為：

式中：ui為速度分量的時均值，ui′為速度分量的脈動值為速度分量乘積的時間平均值。μ為流體的動力粘性系數。

計算中采用的2方程湍流模型為k-ω SST模型。

在數值求解過程中，采用二維非結構三角形網格劃分整個流體域，并采用動網格技術[11]保證物體運動時的網格質量。如圖2所示，為了滿足計算的精確度，網格在半圓柱和水翼周圍分布較密。

圖2 計算網格Fig.2 The mesh grid around semicircular cylinder and foil

圖3 半圓柱振動產生的旋渦Fig.3 The vortices generated by an oscillating semicircular cylinder

3.2 數值計算方法對半圓柱振動問題有效性的驗證

如圖3所示，半圓柱在振動過程中不斷有旋渦脫落到尾流中，隨著旋渦向下游移動，渦強不斷減弱。這些旋渦類似于卡門渦街，在半圓柱尾流場中形成兩排反向旋轉，排列有序的渦列。這些旋渦與卡門渦街的不同之處在于兩排反向旋轉的渦列并不交錯。這兩排渦列在向下游移動的過程中，在與半圓柱直徑平行的方向上總保持一定的距離。計算結果中的現象與Simmons等人的試驗[12]相吻合。在計算過程中，尾流場中的脫落渦在距半圓柱8倍直徑的范圍內可以保持很好的形狀和渦強。擺動水翼總是在這一范圍內運動，因此選擇半圓柱在數值計算中作為旋渦發生裝置是可行的。

采用前面提到的數值方法對半圓柱在不同振幅下的水動力性能進行了計算。為了便于分析計算結果，引入了幾個無因次參數。a=2Ad/d表示半圓柱的無因次振幅，d為半圓柱的直徑。S=s/V0/T、H=h/d分別表示同向旋轉的半圓柱渦之間的無因次距離和反向旋轉的半圓柱渦之間的無因次距離。如圖3所示，其中s為同向旋轉的半圓柱渦之間的實際距離，h為反向旋轉的半圓柱渦之間的實際距離，T為運動周期。

為驗證數值計算方法對于半圓柱振動問題的有效性。對運動參數為f=0.5、V0=2.28 cm/s、d=1.91 cm及a分別取1.0、2.0和3.0這三種工況進行了計算。并與麻省理工學院的實驗結果[13]進行了比較。比較結果如圖4-6所示，計算與實驗符合良好。同時對V0=1.41 cm/s的工況也進行了計算，結果與V0=2.28 cm/s工況近似。

圖4 半圓柱平均阻力比較Fig.4 The comparison of the mean drag force coefficient for an oscillating semicircular cylinder between computational and experimental results

圖5 反向旋轉半圓柱渦無因次距離的比較Fig.5 The comparison of the dimensionless space for vortices in opposite directions between computational and experimental results

圖6 同向旋轉圓柱渦無因次距離的比較Fig.6 The comparison of the dimensionless space for vortices in the same direction between computational and experimental results

圖7 擺動水翼平均推力系數的比較Fig.7 The comparison of the mean thrust force of a flapping foil in the wake of an oscillating semicircular cylinder between computational and experimental results

3.3 數值計算方法對半圓柱水翼串列運動問題有效性的驗證

在水翼運動過程中，首緣和尾緣都存在渦量的集中，從而產生旋渦。當水翼在振動半圓柱尾流中運動時，水翼渦和半圓柱渦相互作用。在計算過程中，通過調整水翼到半圓柱的距離改變水翼渦和半圓柱渦的相互作用模式。圓柱中心和水翼首緣之間的距離表示成無量綱數為σ=l/V0/T，其中l為圓柱中心和水翼首緣之間的實際距離。

為了驗證數值計算方法對半圓柱和水翼串列運動問題的有效性和精確性。對運動參數為f=0.5、V0=1.41cm/s、d=1.91cm、d/C0=0.5、a=2.0、Ad=Af=d、θ0=15°及 σ 分別取 0.75、1.0、1.25 和 1.5 這 4 種工況進行計算。并與麻省理工學院的試驗結果[13]進行了比較。比較結果如圖7所示，計算與試驗結果的變化趨勢基本一致。計算過程中對試驗環境的忽略是造成與試驗結果偏差的主要原因。

4 數值計算結果與分析

本文對運動參數為 f=0.5、V0=1.41 cm/s、d=1.91 cm、d/C0=0.5、a=2.0、Ad=Af=d 及 σ 分別取 0.75、1.0、1.25、1.5和 2.0，水翼搖擺振幅分別取 θ0=15°和 θ0=30°這 10種工況進行了計算。

4.1 旋渦相互作用模式

在水翼上、下表面，半圓柱渦和水翼渦的相互作用模式如表1所示。在計算結果中，半圓柱渦和水翼渦存在4種相互作用模式，這4種旋渦作用模式在其他學者的試驗[14-15]中也有所體現。

表1 旋渦作用模式Tab.1 Modes of action between vortices

模式1：半圓柱渦在向下游移動過程中與水翼渦保持較遠距離，半圓柱渦與水翼渦相互影響較小。如圖8（a）所示，一個旋渦從半圓柱下方脫落。如圖8（b）-8（d）所示，該旋渦在向下游移動過程中，強度不斷減弱，直至消失，在這一過程中半圓柱渦幾乎沒有和水翼渦相互作用。

模式2：半圓柱渦與旋向相同的水翼首緣渦相互作用，并且沿著水翼表面獨立地向下游移動。最終半圓柱渦和水翼首緣渦相繼融入水翼尾緣渦中。如圖9（a）所示，一個半圓柱渦靠近水翼。在圖9（b）中，該半圓柱渦和水翼下方的首緣渦相互作用。如圖9（c）-9（d）所示，半圓柱渦和首緣渦沿水翼表面向水翼尾緣移動，并最終融入水翼尾緣渦中，隨著尾緣渦脫落。

模式3：半圓柱渦在向下游移動的過程中，和水翼首緣渦的相互影響較小，最終融入到相同旋向的水翼尾緣渦中。如圖10（a）所示，一個半圓柱渦靠近水翼，與此同時，位于水翼上方的水翼首緣渦產生，并沿水翼表面向水翼尾緣運動。如圖10（b）所示，半圓柱渦已經到達水翼的首緣，但此時水翼首緣渦已經融入到了水翼尾緣渦中，因此半圓柱渦并沒對水翼首緣渦產生較大影響。如圖10（c）-10（d）中，半圓柱渦融入到水翼尾緣渦中,最終隨尾緣渦脫落。

模式4：半圓柱渦和相反旋向的水翼渦相互作用。如圖11（a）所示，一個半圓柱渦靠近水翼的首緣。在圖11（b）中，半圓柱渦和水翼下方相反旋向的首緣渦相互作用，使首緣渦的強度減弱。如圖11（c）所示，半圓柱渦離開水翼表面，水翼首緣渦融入到水翼尾緣渦中。在圖11（d）中，水翼尾緣渦脫落，半圓柱渦繼續向下游移動。

圖8 模式1Fig.8 Mode 1

圖9 模式2Fig.9 Mode 2

圖10 模式3Fig.10 Mode 3

圖11 模式4Fig.11 Mode 4

4.2 振動半圓柱尾流中的擺動水翼水動力系數

圖12、圖13為10種工況下擺動水翼水動力系數在一個運動周期內隨時間變化的曲線。如文中所述，在這10種工況中存在多種旋渦作用模式，盡管同種水動力系數曲線的形狀相差都不大，但是曲線的幅度變化很大。這說明半圓柱渦對水動力系數隨時間變化規律的影響很小，但對水動力系數幅值的影響很大。在一個周期內，推力系數Cx連續變化兩次，存在兩個峰值，第一個峰值出現在二分之一周期附近，第二個峰值出現在周期結束處。側向力系數Cy在一個周期內連續變化一次，存在兩個大小相等，方向相反的峰值，第一個峰值出現在二分之一周期處，第二個峰值出現在周期結束處。力矩系數Cm的變化規律和側向力系數Cy的變化規律相似，只是與Cy相比，Cm的幅值要小很多。當σ=2.0時半圓柱渦和水翼渦的相互作用很弱（模式1），因此半圓柱渦對水翼水動力系數的影響很小。如圖12所示，在θ0=15°時的水動力系數曲線中，當σ取1.0、1.25和1.5時，水動力曲線的幅值都大于σ=2.0時的幅值。σ=0.75時的水動力曲線幅值小于σ=2.0時的幅值。如圖13所示，在θ0=30°時的水動力系數曲線中，當σ取1.25、1.5時，水動力曲線的幅值都大于σ=2.0時的幅值。當σ取0.75、1.0時，水動力曲線的幅值都小于σ=2.0時的幅值。

圖12 擺動水翼水動力系數（θ0=15°）Fig.12 The hydrodynamic coefficients of a flapping hydrofoil at θ0=15°behind an oscillating semicircular cylinder

圖13 擺動水翼水動力系數（θ0=30°）Fig.13 The hydrodynamic coefficients of a flapping hydrofoil at θ0=30°behind an oscillating semicircular cylinder

4.3 振動半圓柱尾流中的擺動水翼平均推力系數

圖14為振動半圓柱后擺動水翼在θ0=15°、θ0=30°時的平均推力系數曲線。當θ0=15°時，最大平均推力系數Cxm出現在σ=1.25時，比σ=2.0時的Cxm高出約10.8%，半圓柱渦和水翼渦以模式4相互作用。當θ0=30°，最大平均推力系數Cxm=2.02出現在σ=1.25時，比σ=2.0時的Cxm高出約34.7%，半圓柱渦和水翼渦以模式4相互作用。綜上所述，在計算的10種工況中，最大平均推力系數總是發生在σ=1.25時，此時半圓柱渦和相反旋向的水翼渦相互作用（模式4）。

4.4 振動半圓柱尾流中擺動水翼的輸入功率系數和效率

圖15為振動半圓柱后方二維擺動水翼在θ0=15°和θ0=30°時的輸入功率系數曲線和推進效率曲線。當σ=2.0時，半圓柱渦和水翼渦相互作用很小，半圓柱渦對水翼的輸入功率和效率的影響不大。

如圖 15（a）所示，σ=1.25時的輸入功率系數 Pin比σ=2.0時的Pin高出約17.3%，而σ=1.25時的推進效率η卻和σ=2.0時的η接近。這說明水翼在σ=1.25時一部分輸入功率并沒有輸出轉化為推力，而是在半圓柱渦的影響下損失了能量。σ=1.5時的Pin與σ=2.0時的Pin接近，但 σ=1.5時的 η比 σ=2.0時的 η 高出約14.3%。這說明水翼在σ=1.5時除了自身運動提供的輸入功率外，在半圓柱渦的影響下又額外得到了一部分輸入功率，吸收了能量。如圖15（b）所示，σ=1.5時的Pin略小于σ=2.0時的Pin，然而σ=1.5時的η比σ=2.0時的η高出約8.1%，這說明σ=1.5時擺動水翼在半圓柱渦的影響下得到額外的輸入功率，吸收了能量。

圖14 擺動水翼平均推力系數Fig.14 The mean thrust force coefficient Cxmof a flapping hydrofoil

圖15 擺動水翼的輸入功率系數和推進效率Fig.15 The input power coefficient and propulsive efficiency of a flapping hydrofoil

綜上所述，在計算的10種工況中，擺動水翼總是在σ=1.5時吸收能量，此時半圓柱渦和相同旋向的水翼首緣渦相互作用，最終半圓柱渦融入到水翼尾緣渦中（模式2）。

5 結論與展望

本文利用數值方法分析了振動半圓柱尾流中二維擺動水翼的推進性能。計算結果顯示振動半圓柱產生的旋渦對擺動水翼的水動力性能有顯著的影響。可以得到以下結論：

（1）文中的數值計算方法對于分析半圓柱運動和半圓柱、水翼串列運動都是準確、有效的。

（2）在數值計算中，半圓柱可以作為旋渦發生裝置。振動的半圓柱能夠產生兩排旋向相反、排列整齊和高質量的旋渦序列。

（3）半圓柱渦和水翼渦之間存在4種相互作用模式。半圓柱渦對擺動水翼水動力周期變化規律影響較小，對水動力的幅值影響較大。

（4）擺動水翼最大推力系數出現在旋渦相互作用模式4（旋向相同的半圓柱渦和水翼渦相互作用），擺動水翼從半圓柱渦吸收能量的現象出現在旋渦相互作用模式2（旋向相反的半圓柱渦和水翼首緣渦相互作用，半圓柱渦最終融入水翼尾緣渦）。

振動半圓柱產生的旋渦對柔性擺動水翼推進性能的影響將是下一步研究的方向。

[1]George V Lauder,Eliot G Drucker.Morphology and experimental hydrodynamics of fish fin control surfaces[J].IEEE Journal of Oceanic Engineering,2004,29(3):556-571.

[2]Lindsey C C.Form function and locomotory habits in fish.Fish Physiology[M].New York:Academic Press,1978(7):1-100.

[3]Anderson J M,Kerrebrock P A.The Vorticity Control Unmanned Undersea Vehicle(VCUUV):An autonomous robot tuna[C]//Proceedings of the 11th International Symposium on Unmanned Untethered Submersible Technology.Durham,USA,1999:63-70.

[4]Streitlien K,Triantafyllou G S.Efficient foil propulsion through vortex control[J].AIAA Journal,1996,34(11):2315-2318.

[5]Gopalkrishnan R,Triantafyllou M S,Triantafyllou G S,et al.Active vorticity control in a shear flow using a flapping foil[J].Journal of Fluid Mechanics,1994,27(4):1-21.

[6]Beal D N.Propulsion through wake synchronization using flapping-foil[D].Cambridge:Massachusetts Institutes of Technology,2003.

[7]Zhu Q,Wolfgang M J,Yue D K P,Triantafyllou M S.Three-dimensional flow structures and vorticity control in fish-like swimming[J].Journal of Fluid Mech,2002,468:1-28.

[8]Yan Hui,Su Yumin,Yang Liang.Experimentation of fish swimming based on tracking locomotion locus[J].Journal of Bionic Engineering,2008,5(3):258-263.

[9]Triantafyllou M S,Triantafyllou G S,Yue D K P.Hydrodynamics of fishlike swimming[J].Fluid Mech,2000,32:33-53.

[10]成 巍.仿生水下機器人仿真與控制技術研究[D].哈爾濱:哈爾濱工程大學,2004.

[11]王志東,張曉慶,叢文超.柔性擺動水翼弦向變形模式及其對推進性能影響的研究[J].船舶力學,2010,14(7):609-707.Wang Zhidong,Zhang Xiaoqing,Cong Wenchao.Research on the influence of chordwise deflection mode on the propulsion performance of flexible flapping foil[J].Journal of Ship Mechanics,2010,14(7):609-707.

[12]Simmons J E L.Phase-angle measurements between hot-wire signals in the turbulent wake of a two-dimensional bluff body[J].Journal of Fluid Mechanics,1974,64:599-609.

[13]Anderson J M.Vorticity control for efficiency propulsion[D].Cambridge:Massachusetts Institutes of Technology,1996.

[14]Gopalkrishnan R.Vortex induced forces on oscillating bluff cylinders[D].Cambridge:Massachusetts Institutes of Technology,1994.

[15]Streitlien Knut.Extracting energy form unsteady flow through vortex control[D].Cambridge:Massachusetts Institutes of Technology,1994.