艾賢祖， 宋顯成， 趙國(guó)平

(北京精密機(jī)電控制設(shè)備研究所, 北京 100076)

相較于傳統(tǒng)螺旋槳推進(jìn), 魚(yú)類的推進(jìn)方式實(shí)現(xiàn)了槳舵一體, 具有高機(jī)動(dòng)、高效率和低擾動(dòng)等優(yōu)點(diǎn)。而在魚(yú)類眾多的推進(jìn)方式中, 以鲹科魚(yú)類為代表的身體-尾鰭推進(jìn)模式(body and/or caudal fin, BCF)機(jī)動(dòng)性能較好、巡游速度快, 是優(yōu)良的仿生樣本。該模式主要依靠魚(yú)類身體的波動(dòng)及尾部的擺動(dòng)產(chǎn)生推進(jìn)力, 其中鲹科魚(yú)類通常通過(guò)短時(shí)爆發(fā)、 長(zhǎng)時(shí)滑行來(lái)減少效率損失[1], 在這個(gè)過(guò)程中尾部的擺動(dòng)將直接影響推進(jìn)性能的好壞。尤其是作為魚(yú)尾主要部分的尾鰭, 其運(yùn)動(dòng)參數(shù)對(duì)推進(jìn)性能有著重要影響。研究表明[2], 相較于升沉幅值,搖擺幅度對(duì)尾鰭推進(jìn)性能的影響要大得多。因此研究擺幅對(duì)尾鰭推進(jìn)性能的影響有著更重要意義。

針對(duì)尾鰭運(yùn)動(dòng)參數(shù)對(duì)其推進(jìn)性能影響,許多學(xué)者開(kāi)展了研究。俞經(jīng)虎等[3]通過(guò)一種三節(jié)仿生機(jī)器魚(yú)模型,分別計(jì)算水動(dòng)力和慣性力引起的尾鰭驅(qū)動(dòng)力矩的值。劉軍考等[4]通過(guò)實(shí)驗(yàn)研究表明, 在一定擺動(dòng)頻率范圍內(nèi), 推進(jìn)速度與擺動(dòng)頻率在一定范圍內(nèi)近似成正比關(guān)系, 而當(dāng)擺動(dòng)頻率一定時(shí), 推進(jìn)速度與擺幅在一定范圍內(nèi)近似成正比關(guān)系。張曦等[5]通過(guò)數(shù)值計(jì)算及實(shí)驗(yàn)得出, 隨著斯特勞哈爾數(shù)增大,仿金槍魚(yú)尾鰭的平均推力系數(shù)增大,推進(jìn)效率降低。劉葳興等[6]以能量的視角提出一種對(duì)剛性和柔性拍動(dòng)翼都適用的推進(jìn)效率數(shù)值計(jì)算方法, 并系統(tǒng)地研究了各參數(shù)對(duì)三維柔性翼推進(jìn)性能的影響。胡健等[7]采用計(jì)算流體力學(xué)方法分析了擺動(dòng)水翼的水動(dòng)力特性, 詳細(xì)分析了翼型的水動(dòng)力性能隨斯特勞哈爾數(shù)(Strouhal number)Sr的變化的影響。

對(duì)比其他形狀尾鰭, 金槍魚(yú)的月牙形尾鰭具有極其出色的推進(jìn)性能[8],有著可預(yù)見(jiàn)的廣闊應(yīng)用前景。同時(shí),鑒于擺幅對(duì)其推進(jìn)性能影響的重要性,以金槍魚(yú)尾鰭形狀特征為分析對(duì)象, 基于浸入邊界法(immersed boundary method,IBM)開(kāi)展擺幅對(duì)其推進(jìn)性能影響研究。首先,對(duì)金槍魚(yú)尾鰭外形特征以及運(yùn)動(dòng)特征進(jìn)行數(shù)學(xué)描述,并建立模型。然后,以蘇玉民等[2]進(jìn)行的擺動(dòng)尾鰭試驗(yàn)為算例,對(duì)比驗(yàn)證了浸入邊界數(shù)值方法的可行性。最后,在給定斯特勞哈爾數(shù)條件下, 對(duì)比分析仿金槍魚(yú)尾鰭在不同搖擺幅度下推進(jìn)性能的變化,并提煉相關(guān)規(guī)律，為水下仿金槍魚(yú)尾鰭推進(jìn)機(jī)構(gòu)的設(shè)計(jì)和控制提供參考。

1 數(shù)值方法驗(yàn)證

1.1 浸入邊界法簡(jiǎn)介

浸入邊界法(IBM) 是采用歐拉變量去描述流體的動(dòng)態(tài)、利用拉格朗日變量描述結(jié)構(gòu)的運(yùn)動(dòng)邊界、用光滑Delta近似函數(shù)通過(guò)分布節(jié)點(diǎn)力和插值速度來(lái)表示流場(chǎng)和結(jié)構(gòu)物的交互作用的方法[9]。它整個(gè)流場(chǎng)計(jì)算都使用笛卡爾網(wǎng)格, 無(wú)需處理從物理平面到計(jì)算平面的坐標(biāo)和網(wǎng)格轉(zhuǎn)換問(wèn)題, 因而可以大大提高計(jì)算效率，而且節(jié)省了網(wǎng)格生成所需的時(shí)間[10]。

浸入實(shí)體 (immersed solid)功能是以浸入邊界法為原理, 允許用戶對(duì)一些涉及剛體對(duì)象穿過(guò)流體域的運(yùn)動(dòng)進(jìn)行穩(wěn)態(tài)或瞬態(tài)模擬。利用該功能, 可以不用進(jìn)行幾何重構(gòu)和網(wǎng)格重新劃分而獲得將單個(gè)連續(xù)流體域分割成多個(gè)不連續(xù)流體域的流體仿真結(jié)果。而在以前的流體分析中, 一般要求計(jì)算的流體域是連續(xù)的, 即流體計(jì)算需要計(jì)算連續(xù)性方程。浸入實(shí)體功能可以以小幅犧牲精度為代價(jià), 大大加快計(jì)算速度及計(jì)算穩(wěn)定性。

1.2 方法驗(yàn)證

首先,需驗(yàn)證IBM數(shù)值方法的有效性。選取蘇玉民等[2]在循環(huán)水槽中所進(jìn)行的擺動(dòng)尾鰭水動(dòng)力性能研究試驗(yàn)作為對(duì)比對(duì)象, 對(duì)其模型進(jìn)行數(shù)值仿真,最后比對(duì)數(shù)值結(jié)果與試驗(yàn)數(shù)據(jù)來(lái)驗(yàn)證方法。仿真模型如圖1所示。模型特征弦長(zhǎng)C0為0.3 m。針對(duì)外流場(chǎng),通常認(rèn)為上下游跨距超過(guò)10倍特征長(zhǎng)度即可認(rèn)為遠(yuǎn)場(chǎng)無(wú)擾動(dòng),側(cè)向跨距超過(guò)5倍特征長(zhǎng)度即可認(rèn)為遠(yuǎn)場(chǎng)無(wú)擾動(dòng)。因此,仿真的計(jì)算域取x軸方向(上下游方向)10倍弦長(zhǎng),在y軸方向取6倍弦長(zhǎng),在z軸方向取6倍弦長(zhǎng),外邊界設(shè)置為壓力遠(yuǎn)場(chǎng)邊界,入口給定流速入口,出口給定壓力出口邊界,采用kω-SST湍流模型,最終的網(wǎng)格無(wú)關(guān)性網(wǎng)格數(shù)為300萬(wàn)。

圖1 方法驗(yàn)證數(shù)值模擬模型Fig.1 Method validation numerical simulation model

通過(guò)ANSYS 流體動(dòng)力學(xué)仿真軟件(ANSYS CFX, CFX)中浸入實(shí)體功能對(duì)擺動(dòng)尾鰭試驗(yàn)?zāi)Ｐ瓦M(jìn)行數(shù)值計(jì)算,得到一個(gè)周期內(nèi)的推力系數(shù)Cx。將試驗(yàn)結(jié)果與數(shù)值模擬結(jié)果進(jìn)行對(duì)比, 如圖2所示?？梢钥闯?模擬結(jié)果與試驗(yàn)數(shù)據(jù)曲線變化趨勢(shì)大體一致,只在局部的波谷處有較大偏差。這是由于浸入邊界法本質(zhì)上是基于固壁邊界附面層采用源項(xiàng)等效的一種近似,這種處理在層流及小分離流動(dòng)時(shí)能達(dá)到較好的模擬精度,但當(dāng)局部出現(xiàn)瞬時(shí)大分離流動(dòng)時(shí),由于湍流固有的結(jié)構(gòu)復(fù)雜性和不確定性的加劇,此時(shí)單純的源項(xiàng)是不足以準(zhǔn)確近似湍流的,必然會(huì)出現(xiàn)較大數(shù)值偏差,但受主流區(qū)流動(dòng)的制約,這種局部偏差有限,不會(huì)對(duì)流場(chǎng)發(fā)展趨勢(shì)有質(zhì)的改變,以其作為推進(jìn)性能的定性判斷還是可以作為參考的。因此,IBM方法對(duì)本文問(wèn)題研究具有有效性。

圖2 計(jì)算結(jié)果與算例中試驗(yàn)結(jié)果比較Fig.2 Comparison of calculation results with test results in calculation examples

2 仿金槍魚(yú)尾鰭幾何建模

在新月形尾鰭外形參數(shù)方面, 根據(jù)陳宏[11]在設(shè)計(jì)一種金槍魚(yú)外形的仿生機(jī)器魚(yú)過(guò)程中, 通過(guò)最小二乘法對(duì)真實(shí)金槍魚(yú)輪廓進(jìn)行擬合得到的曲線公式對(duì)新月形尾鰭進(jìn)行設(shè)計(jì)。根據(jù)其給出的公式, 可以得到本文中的新月形尾鰭曲線公式為

式(1)中:x為魚(yú)體沿著體長(zhǎng)方向坐標(biāo)；y為魚(yú)體沿著體高方向坐標(biāo)；lp為尾柄最細(xì)位置到魚(yú)頭頂點(diǎn)位置的長(zhǎng)度, 為1 700 mm；lm為魚(yú)體在中心線上總長(zhǎng)度, 為1 828 mm；l為魚(yú)體全長(zhǎng), 為2 000 mm。

根據(jù)式(1)中給出的真實(shí)金槍魚(yú)魚(yú)體輪廓曲線函數(shù), 可以得到金槍魚(yú)垂直基準(zhǔn)面外形輪廓線, 如圖3所示。

圖3 金槍魚(yú)垂直基準(zhǔn)面輪廓線Fig.3 Vertical plane contour of tuna

根據(jù)上述金槍魚(yú)體輪廓曲線函數(shù), 對(duì)新月形尾鰭外形參數(shù)進(jìn)行近似選擇, 圖4(a)所示為新月形尾鰭外形參數(shù)，其三維建模如圖4(b)所示。根據(jù)式(1)進(jìn)行計(jì)算, 選擇并優(yōu)化新月形尾鰭外形參數(shù), 可以給出尾鰭長(zhǎng)度lt為300 mm, 鰭梢長(zhǎng)度lf為80 mm, 中弦長(zhǎng)lc為128 mm, 后掠角α=47.4°, 展長(zhǎng)b為242 mm。

在三維建模過(guò)程中, 根據(jù)新月形尾鰭形狀及其主要參數(shù)建立模型, 同時(shí)新月形尾鰭的剖面翼型選擇NACA0015系列翼型, 進(jìn)行三維建?？梢缘玫叫略滦挝馋?。

圖4 新月外形參數(shù)及三維模型Fig.4 Lunate shape parameters and 3D model

3 仿金槍魚(yú)尾鰭不同擺幅流場(chǎng)模擬

3.1 擺動(dòng)尾鰭運(yùn)動(dòng)參數(shù)定義

對(duì)于以金槍魚(yú)為代表的身體-尾鰭推進(jìn)模式(body and/or caudal fin, BCF), 其主要推進(jìn)力來(lái)自于后1/3部分。參照?qǐng)D5可以將尾鰭進(jìn)行的復(fù)雜運(yùn)動(dòng)簡(jiǎn)化為以下3個(gè)簡(jiǎn)單運(yùn)動(dòng)。

圖5 金槍魚(yú)尾鰭運(yùn)動(dòng)軌跡Fig.5 Movement track of a tuna-like caudal fin

(1)魚(yú)在向前游動(dòng)過(guò)程中, 尾鰭在X軸方向上及前進(jìn)方向上產(chǎn)生的平動(dòng)運(yùn)動(dòng), 即巡航速度V0。

(2)尾鰭跟隨魚(yú)尾擺動(dòng)進(jìn)行的平動(dòng)運(yùn)動(dòng), 即升沉運(yùn)動(dòng), 沿Y軸方向可表示為

y(t)=Asin(2πft) (2)

(3)尾鰭繞自身擺動(dòng)軸的擺動(dòng)運(yùn)動(dòng), 即搖擺運(yùn)動(dòng), 繞Z軸擺動(dòng)的角位移可表示為

θ(t)=θ0sin(2πft-φ) (3)

式(2)、式(3)中:A為尾鰭的升沉運(yùn)動(dòng)的幅值,m；f為耦合運(yùn)動(dòng)的頻率,Hz；θ0為尾鰭的搖擺運(yùn)動(dòng)幅度；φ為升沉運(yùn)動(dòng)和搖擺運(yùn)動(dòng)相位差。

對(duì)式(2)、式(3)進(jìn)行求導(dǎo), 可以得到尾鰭沿Y軸方向升沉運(yùn)動(dòng)的速度和繞Z軸搖擺運(yùn)動(dòng)的角速度分別為

V(t)=2πfAcos(2πft) (4)

ω(t)=2πfθ0cos(2πft-φ) (5)

3.2 水動(dòng)力及相關(guān)系數(shù)定義

(1)對(duì)推進(jìn)力進(jìn)行無(wú)量綱化處理得到推力系數(shù)Cx,可以表示為

式(6)中:Fx為推進(jìn)力,N；ρ為水的密度, kg/m3；S為尾鰭投影面積,m2。

(2) 斯特勞哈爾數(shù)Sr是影響擺動(dòng)尾鰭推進(jìn)性能重要的無(wú)量綱參數(shù),表示為

Sr=fb/|V0| (7)

式(7)中:b為尾鰭運(yùn)動(dòng)過(guò)程中脫落尾渦的寬度, 近似表達(dá)為升沉運(yùn)動(dòng)幅值的2倍。

(3)擺動(dòng)尾鰭的推進(jìn)效率可以表示為

式(8)中:Cxm為平均推力系數(shù),意義為Cx在一個(gè)周期內(nèi)的平均值；Cy為側(cè)向力系數(shù)；Cm為力矩系數(shù)。

3.3 變擺幅的數(shù)值模擬實(shí)現(xiàn)

根據(jù)計(jì)算數(shù)據(jù)進(jìn)行三維建模, 分別用六面體、四面體對(duì)流域、尾鰭模型進(jìn)行網(wǎng)格劃分,計(jì)算模型與圖1相似。通過(guò)浸入實(shí)體功能對(duì)剛體模型尾鰭穿過(guò)以水為介質(zhì)的流體域的擺動(dòng)過(guò)程進(jìn)行模擬, 從而得到尾鰭在恒流, 即魚(yú)在巡航狀態(tài)下, 受到的推進(jìn)力等相關(guān)參數(shù)。

改變尾鰭的搖擺運(yùn)動(dòng)幅度θ0, 使其依次為20°、25°、30°、35°,在這4種情況下對(duì)比分析搖擺運(yùn)動(dòng)幅度對(duì)尾鰭推進(jìn)性能的影響。其余運(yùn)動(dòng)參數(shù)設(shè)置為周期T為1.6 s, 相位差φ為-50°, 尾鰭的升沉運(yùn)動(dòng)的幅值A(chǔ)為0.2 m,V0為0.3 m/s, 對(duì)應(yīng)的Sr為0.4。

設(shè)置剛體尾鰭類型為浸入實(shí)體,選取計(jì)算時(shí)間為8 s, 時(shí)間步長(zhǎng)為0.01 s, 選擇剛體尾鰭運(yùn)動(dòng)方式為規(guī)定位移, 將其運(yùn)動(dòng)分解為平移與旋轉(zhuǎn)進(jìn)行設(shè)置。

4 結(jié)果與分析

基于浸入實(shí)體功能對(duì)尾鰭擺動(dòng)運(yùn)動(dòng)進(jìn)行仿真模擬, 得到壓力分布、尾渦分布以及各個(gè)時(shí)刻的推進(jìn)力大小, 并進(jìn)一步獲取推力系數(shù)以及效率。

4.1 尾鰭表面壓力分布

圖6所示是在θ0為20°的情況下, 尾鰭表面的壓力分布云圖。從圖6中可以看出，尾鰭正面的壓力分布與1/2周期后的背面壓力分布基本相同,尾鰭正反兩面從高壓變?yōu)榈蛪?再?gòu)牡蛪恨D(zhuǎn)變?yōu)楦邏? 相互對(duì)應(yīng)著按周期規(guī)律變化,由此可以看出尾鰭表面的壓力分布是呈周期性規(guī)律進(jìn)行變化的。

圖6 尾鰭表面壓力分布在一個(gè)周期內(nèi)的變化Fig.6 Changes in pressure distribution on the surface of the caudal fin over a period

4.2 尾渦分布

尾鰭從最大振幅處向平衡位置擺動(dòng)過(guò)程中, 由于尾鰭擺動(dòng)形成的低壓區(qū)會(huì)逐漸向尾鰭鰭尖處移動(dòng), 并最終脫落形成尾渦。這些從尾鰭上脫落下來(lái)的尾渦, 以反卡門渦街的形式進(jìn)行排列, 使尾鰭后面的流場(chǎng)中流體呈噴射流動(dòng)狀態(tài)。在這些噴射流體的反作用下, 尾鰭獲得向前的推進(jìn)力。因此, 在尾鰭擺動(dòng)產(chǎn)生推進(jìn)力的過(guò)程中尾渦有著不可替代的重要作用。在三維流場(chǎng)中, 尾渦的形狀和分布非常復(fù)雜, 通過(guò)ANSYS 流體動(dòng)力學(xué)仿真后處理軟件(ANSYS CFD-POST, CFD-POST)中提供的對(duì)物理量進(jìn)行三維體積渲染可以直觀地顯示流場(chǎng)中的尾渦以及渦量強(qiáng)度。

尾渦渦量的體積渲染圖圖形遮蓋嚴(yán)重, 這里取Z=0的平面進(jìn)行對(duì)尾渦的觀察分析。通過(guò)圖7可以看出尾渦的脫落規(guī)律及渦量大小的變化。尾渦在尾鰭由最大振幅處向平衡位置擺動(dòng)過(guò)程中形成, 故在一個(gè)運(yùn)動(dòng)周期內(nèi)會(huì)脫落生成兩次, 并在脫落后不斷擴(kuò)散, 渦量不斷減小, 一周期內(nèi)先后脫落的尾渦在尾鰭后按次序分布排列, 最終所有尾渦以反卡門渦街形式排列。

圖7 一個(gè)周期內(nèi)尾渦的變化Fig.7 Variation of wake vortex in a cycle

4.3 推進(jìn)系數(shù)及推進(jìn)性能影響

通過(guò)浸入實(shí)體功能可以計(jì)算得到仿金槍魚(yú)尾鰭受到流體的作用力。數(shù)值計(jì)算后得到每個(gè)時(shí)間步下在尾鰭受到的推進(jìn)力, 從而計(jì)算得到推進(jìn)力系數(shù)變化如圖8所示。推進(jìn)系數(shù)Cx對(duì)擺動(dòng)周期比t/T按照類正弦規(guī)律周期性變化, 隨著搖擺幅度θ0增大, 其幅值也在不斷變大, 因此增大搖擺幅度θ0, 使尾鰭與流體間作用力變大。

圖8 一個(gè)運(yùn)動(dòng)周期內(nèi)不同搖擺幅度下推進(jìn)系數(shù)的變化Fig.8 Variation of propulsion coefficient under different swing amplitudes during a motion cycle

進(jìn)一步對(duì)數(shù)據(jù)進(jìn)行處理, 可以得到平均推力系數(shù)Cxm及效率η隨搖擺幅度θ0變化曲線,如圖9所示。平均推力系數(shù)Cxm隨搖擺幅度θ0變大不斷變大, 同時(shí)效率η隨著搖擺幅度θ0變大不斷變小。這表明隨著搖擺幅度θ0的增加, 尾鰭產(chǎn)生的推進(jìn)力變大, 但其推進(jìn)效率不斷下降, 這也符合魚(yú)類啟動(dòng)過(guò)程中擺動(dòng)幅度大、獲得加速度大、但效率遠(yuǎn)低于巡航狀態(tài)的特點(diǎn)。

圖9 搖擺幅度對(duì)尾鰭推進(jìn)性能的影響Fig.9 Effect of swing amplitude on the tail fin propulsion performance

5 結(jié)論

基于浸入邊界法研究了仿金槍魚(yú)尾鰭不同搖擺幅度θ0運(yùn)動(dòng)下的推進(jìn)性能, 得到符合實(shí)際的結(jié)果。

(1)通過(guò)仿金槍魚(yú)尾鰭運(yùn)動(dòng)一個(gè)周期內(nèi)壓力場(chǎng)云圖、尾渦變化圖, 分析了在擺動(dòng)過(guò)程中尾鰭兩側(cè)壓力的變化規(guī)律以及尾渦的變化規(guī)律;在Sr恒定的流場(chǎng)中, 推力系數(shù)曲線周期性變化, 呈類正弦曲線形狀;隨著搖擺幅度θ0的增加, 尾鰭產(chǎn)生的推進(jìn)力變大, 但其推進(jìn)效率不斷下降。

(2)研究驗(yàn)證了IBM方法模擬尾翼非定常流動(dòng)的有效性, 實(shí)現(xiàn)了仿金槍魚(yú)尾鰭在流場(chǎng)中的復(fù)合運(yùn)動(dòng)進(jìn)行模擬, 獲取了尾鰭表面壓力分布及尾渦渦量的分布,得到不同尾鰭搖擺幅度對(duì)推力系數(shù)及效率的影響差別, 提取了利于提高推進(jìn)性能的擺幅影響規(guī)律,對(duì)仿金槍魚(yú)尾鰭機(jī)構(gòu)的設(shè)計(jì)和控制具有重要參考意義。