張 曦，蘇玉民，王兆立

（哈爾濱工程大學(xué) 水下智能機(jī)器人國防科技重點(diǎn)實(shí)驗(yàn)室，哈爾濱 150001）

1 引 言

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

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

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

2 數(shù)值計(jì)算模型

如圖1所示，半圓柱通過沿y軸升沉運(yùn)動(dòng)產(chǎn)生旋渦。半圓柱振動(dòng)滿足的運(yùn)動(dòng)規(guī)律可以表示為：

式中：Ad為半圓柱升沉運(yùn)動(dòng)的振幅，f為半圓柱的運(yùn)動(dòng)頻率。

如圖1所示，二維擺動(dòng)水翼既沿著y軸升沉又繞著水翼首緣點(diǎn)搖擺。因此水翼的擺動(dòng)運(yùn)動(dòng)規(guī)律[8]可以表示為：

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

圖1 半圓柱和水翼的周期性運(yùn)動(dòng)Fig.1 The cycle motion of an oscillating semicircular and a flapping hydrofoil

擺動(dòng)水翼的推力系數(shù)Cx，側(cè)向力系數(shù)Cy和以水翼首緣點(diǎn)為作用點(diǎn)的力矩系數(shù)Cm由相應(yīng)的推力Cx、側(cè)向力Cy和力矩Cm沿水翼表面積分[9]確定：

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

擺動(dòng)水翼的推進(jìn)效率[10]表示為：

式中：Cxm為推力系數(shù)Cx在一個(gè)周期內(nèi)的平均值。

輸入功率系數(shù)表示為：

3 數(shù)值計(jì)算方法

3.1 控制方程

對于不可壓縮的粘性流動(dòng)，在笛卡爾坐標(biāo)系下，用張量的形式表示的時(shí)均連續(xù)性方程和RANS方程可以寫為：

式中：ui為速度分量的時(shí)均值，ui′為速度分量的脈動(dòng)值為速度分量乘積的時(shí)間平均值。μ為流體的動(dòng)力粘性系數(shù)。

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

在數(shù)值求解過程中，采用二維非結(jié)構(gòu)三角形網(wǎng)格劃分整個(gè)流體域，并采用動(dòng)網(wǎng)格技術(shù)[11]保證物體運(yùn)動(dòng)時(shí)的網(wǎng)格質(zhì)量。如圖2所示，為了滿足計(jì)算的精確度，網(wǎng)格在半圓柱和水翼周圍分布較密。

圖2 計(jì)算網(wǎng)格Fig.2 The mesh grid around semicircular cylinder and foil

圖3 半圓柱振動(dòng)產(chǎn)生的旋渦Fig.3 The vortices generated by an oscillating semicircular cylinder

3.2 數(shù)值計(jì)算方法對半圓柱振動(dòng)問題有效性的驗(yàn)證

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

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

為驗(yàn)證數(shù)值計(jì)算方法對于半圓柱振動(dòng)問題的有效性。對運(yùn)動(dòng)參數(shù)為f=0.5、V0=2.28 cm/s、d=1.91 cm及a分別取1.0、2.0和3.0這三種工況進(jìn)行了計(jì)算。并與麻省理工學(xué)院的實(shí)驗(yàn)結(jié)果[13]進(jìn)行了比較。比較結(jié)果如圖4-6所示，計(jì)算與實(shí)驗(yàn)符合良好。同時(shí)對V0=1.41 cm/s的工況也進(jìn)行了計(jì)算，結(jié)果與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 反向旋轉(zhuǎn)半圓柱渦無因次距離的比較Fig.5 The comparison of the dimensionless space for vortices in opposite directions between computational and experimental results

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

圖7 擺動(dòng)水翼平均推力系數(shù)的比較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 數(shù)值計(jì)算方法對半圓柱水翼串列運(yùn)動(dòng)問題有效性的驗(yàn)證

在水翼運(yùn)動(dòng)過程中，首緣和尾緣都存在渦量的集中，從而產(chǎn)生旋渦。當(dāng)水翼在振動(dòng)半圓柱尾流中運(yùn)動(dòng)時(shí)，水翼渦和半圓柱渦相互作用。在計(jì)算過程中，通過調(diào)整水翼到半圓柱的距離改變水翼渦和半圓柱渦的相互作用模式。圓柱中心和水翼首緣之間的距離表示成無量綱數(shù)為σ=l/V0/T，其中l(wèi)為圓柱中心和水翼首緣之間的實(shí)際距離。

為了驗(yàn)證數(shù)值計(jì)算方法對半圓柱和水翼串列運(yùn)動(dòng)問題的有效性和精確性。對運(yùn)動(dòng)參數(shù)為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 種工況進(jìn)行計(jì)算。并與麻省理工學(xué)院的試驗(yàn)結(jié)果[13]進(jìn)行了比較。比較結(jié)果如圖7所示，計(jì)算與試驗(yàn)結(jié)果的變化趨勢基本一致。計(jì)算過程中對試驗(yàn)環(huán)境的忽略是造成與試驗(yàn)結(jié)果偏差的主要原因。

4 數(shù)值計(jì)算結(jié)果與分析

本文對運(yùn)動(dòng)參數(shù)為 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種工況進(jìn)行了計(jì)算。

4.1 旋渦相互作用模式

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

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

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

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

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

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

圖8 模式1Fig.8 Mode 1

圖9 模式2Fig.9 Mode 2

圖10 模式3Fig.10 Mode 3

圖11 模式4Fig.11 Mode 4

4.2 振動(dòng)半圓柱尾流中的擺動(dòng)水翼水動(dòng)力系數(shù)

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

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

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

4.3 振動(dòng)半圓柱尾流中的擺動(dòng)水翼平均推力系數(shù)

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

4.4 振動(dòng)半圓柱尾流中擺動(dòng)水翼的輸入功率系數(shù)和效率

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

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

圖14 擺動(dòng)水翼平均推力系數(shù)Fig.14 The mean thrust force coefficient Cxmof a flapping hydrofoil

圖15 擺動(dòng)水翼的輸入功率系數(shù)和推進(jìn)效率Fig.15 The input power coefficient and propulsive efficiency of a flapping hydrofoil

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

5 結(jié)論與展望

本文利用數(shù)值方法分析了振動(dòng)半圓柱尾流中二維擺動(dòng)水翼的推進(jìn)性能。計(jì)算結(jié)果顯示振動(dòng)半圓柱產(chǎn)生的旋渦對擺動(dòng)水翼的水動(dòng)力性能有顯著的影響。可以得到以下結(jié)論：

（1）文中的數(shù)值計(jì)算方法對于分析半圓柱運(yùn)動(dòng)和半圓柱、水翼串列運(yùn)動(dòng)都是準(zhǔn)確、有效的。

（2）在數(shù)值計(jì)算中，半圓柱可以作為旋渦發(fā)生裝置。振動(dòng)的半圓柱能夠產(chǎn)生兩排旋向相反、排列整齊和高質(zhì)量的旋渦序列。

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

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

振動(dòng)半圓柱產(chǎn)生的旋渦對柔性擺動(dòng)水翼推進(jìn)性能的影響將是下一步研究的方向。

[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]成 巍.仿生水下機(jī)器人仿真與控制技術(shù)研究[D].哈爾濱:哈爾濱工程大學(xué),2004.

[11]王志東,張曉慶,叢文超.柔性擺動(dòng)水翼弦向變形模式及其對推進(jìn)性能影響的研究[J].船舶力學(xué),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.