李寧宇,蘇廣勝,趙云鶴,孫向東,張嗣祺,蘇玉民
(1. 哈爾濱工程大學 水下機器人技術(shù)重點實驗室,黑龍江 哈爾濱 150001; 2. 中國石油天然氣管道工程有限公司,河北 廊坊 065000; 3. 杭州和利時自動化有限公司,浙江 杭州 310000)
近年來,利用魚類游動的推進機理所研制的仿生水下機器人在海底地形掃描、海洋資源開發(fā)等領(lǐng)域的應(yīng)用日益增多。就游動時所依靠的生理部位來講,可分為身體和/或尾鰭(body and/or caudal fin, 簡稱BCF)游動模式以及中央鰭和/或?qū)?median and/or paired fin, 簡稱MPF)游動模式[1]。由于相對簡單的設(shè)計和較高的推進性能,當前的仿生推進系統(tǒng)多采用BCF模式[2-3],然而不少種類的魚采用MPF模式進行姿態(tài)控制、提供輔助推進力等,以及在低速前進時直接將胸鰭作為主推進器[4],海龜、企鵝等也使用類胸鰭的鰭狀肢來推進身體向前游動。
有關(guān)仿生拍動翼(可看作胸鰭的一個簡化模型)的研究發(fā)現(xiàn),反卡門渦街為二維拍動翼尾流場的重要特征[5-6],然而對于三維翼而言,其尾流場則由兩個傾斜渦列所支配[7-8]。而胸鰭的形狀和運動比三維翼更加復(fù)雜,胸鰭尾渦結(jié)構(gòu)的關(guān)鍵特征究竟是什么還有待研究。
Suzuki等[9]開發(fā)了仿黑鱸魚機械胸鰭推進系統(tǒng),并開展實驗研究工作,分析胸鰭的瞬時以及時均水動力特性。Li等[10]利用FLUENT軟件及其動網(wǎng)格技術(shù)分析了胸鰭非定常運動尾流場中流向、展向系列切片上的渦的脫落和耗散過程。Lauder等[11]和Bozkurttas等[12]試驗觀測了太陽魚穩(wěn)定游動狀態(tài)下的胸鰭幾何和運動學數(shù)據(jù),并進行與試驗相匹配的數(shù)值模擬以提供詳細的流場和水動力定量信息。陳宏等[13]建立了在胸鰭擺動模式下仿生機器人游動的物理模型,利用理論方法計算巡游時所受流體力,并進行數(shù)值仿真分析。Ramamurti等[14]對瀨魚胸鰭的游動行為進行數(shù)字化,并實施了相應(yīng)的三維準定常和非定常胸鰭繞流數(shù)值研究。蘇玉民等[15]開發(fā)了非定常渦格法程序,用行程時間比來模型化胸鰭的運動,計算分析了仿生胸鰭的推進性能。Hu[16]研究了在基于阻力的游動模式下胸鰭的水動力學問題,利用CFD計算給出了胸鰭的推力、推進效率和二維流場渦量分布。上述研究在仿生胸鰭的水動力性能預(yù)報和流場定量信息的獲取方面取得了一些有價值的成果,但胸鰭三維尾渦結(jié)構(gòu)的非定常演化過程還沒有被深入地探討,并且缺乏對胸鰭運動參數(shù)影響其水動力性能的機理的理解。
首先參考活體魚胸鰭試驗和運動學分析的有關(guān)文獻進行仿生胸鰭計算模型的建立,然后通過與模型試驗和其他數(shù)值方法的比較,驗證所采用的數(shù)值方法模擬胸鰭繞流問題的精確性和可靠性,接下來研究仿生鰭瞬態(tài)運動過程中的尾渦結(jié)構(gòu)演化與非定常水動力的關(guān)系,并探明支配胸鰭尾流場的關(guān)鍵渦結(jié)構(gòu),最后關(guān)注運動參數(shù)對仿生鰭的水動力性能和流場渦結(jié)構(gòu)的影響,并根據(jù)研究成果對仿胸鰭推進系統(tǒng)的設(shè)計提出合理化建議。
參考真實胸鰭的游動試驗[17-18],CFD計算中的仿生鰭幾何模型如圖1所示,其中鰭的最大弦長C為0.155 m,展弦比約為1.2,仿生鰭的尺寸約為試驗中黑鱸魚的胸鰭的6.2倍。
圖1 仿生胸鰭的模型Fig. 1 Model of a bio-inspired pectoral fin
隆頭魚科模式魚類胸鰭的典型游動模式有基于升力的游動模式和基于阻力的游動模式,其中前者由縱傾和上下拍動運動組成,而后者由縱傾和前后拍動運動組成[4]以及Kato[18]對黑鱸魚胸鰭的觀測和運動學分析,這種魚的胸鰭可作為基于阻力的游動模式的典型。胸鰭的前后拍動運動可以表達為
φR(t)=φRC-φRAcos(2πf)
(1)
式中:φRA為前后拍動運動的幅值,φRC為前后拍動運動的平均值,f為運動頻率,t為時間。
胸鰭的縱傾運動可由下式給出:
φFE(t)=φFEC-φFEAcos(2πf+ΔφFE)
(2)
式中:φFEA為縱傾運動的幅值,φFEC為縱傾運動的平均值,ΔφFE表示前后拍動與縱傾運動之間的相位差。圖2中給出了在一個周期內(nèi)胸鰭非定常運動的8個典型瞬時,為便于表明在胸鰭運動過程中鰭和魚體的相對位置的變化,魚體也被繪制出來以更清楚地顯示胸鰭的運動情況(魚體僅為顯示目的并不包含在當前的數(shù)值計算中)。
圖2 在一個周期內(nèi)的胸鰭運動情況Fig. 2 Motion of the pectoral fin during one period
仿生拍動翼/鰭的研究中常用的無量綱參數(shù),雷諾數(shù)(Re)和斯特勞哈爾數(shù)(St),定義如下:
Re=UC/ν
(3)
式中:ν為流體的運動黏性系數(shù)。
St=fC/U
(4)
基于所建立的坐標系統(tǒng),胸鰭的推力系數(shù)CT和升力系數(shù)CL可以表示為
(5)
式中:ρ為流體密度,F(xiàn)x為推力,-Fz表示升力(考慮到z軸向下),Aplan為胸鰭的投影面積。對CT和CL分別在一個運動周期內(nèi)取平均,有:
(6)
式中:CTA為平均推力系數(shù),CLA為平均升力系數(shù),T表示周期。
推進效率為胸鰭的平均輸出功率Pout與平均輸入功率Pin的比值,可用來流速度與平均推力的乘積來計算Pout,而Pin通過水對胸鰭做負功而被水吸收,最終得到計算胸鰭的推進效率的公式為
(7)
鑒于在本研究中所應(yīng)用的數(shù)值方法(一種浸入式邊界方法)已在文獻[19-21]中被詳細闡述,這里僅對該方法的一些基本特征加以簡要描述??刂菩伥捓@流問題的方程為非定常、不可壓縮、三維Navier-Stokes方程:
(8)
式中:u和p分別表示流體的速度和壓力,ρ為流體的密度。采用有限體積法離散。一個分裂步方法(fractional-step method)[22]被用于壓力速度耦合方程組的求解,流場網(wǎng)格劃分形式為規(guī)則的六面體網(wǎng)格,物面網(wǎng)格劃分形式為三角形網(wǎng)格。浸沒邊界對流體的影響通過引入虛擬單元(ghost-cell)來表示,這種方法的主要思想是利用流場的局部插值重構(gòu)來滿足物面的壓力和速度邊界條件。
在文獻[19-21]中,已通過若干算例對上述數(shù)值方法進行了廣泛的驗證,包括繞固定球體的流動、做大幅度運動的拍動翼和均勻流中圓柱體受迫振蕩。本節(jié)針對所關(guān)注的仿生胸鰭水動力性能數(shù)值模擬,首先進行網(wǎng)格和計算域的無關(guān)性研究,然后將浸入邊界法計算結(jié)果與試驗結(jié)果和前人采用其他數(shù)值方法得到的結(jié)果進行比較,進一步驗證本文數(shù)值方法的正確性。
數(shù)值計算中所采用的標準網(wǎng)格如圖3所示。胸鰭表面由13 732個三角形網(wǎng)格離散;背景流場由693萬個笛卡爾網(wǎng)格離散,為清晰顯示流場中的網(wǎng)格分布細節(jié),三個相互垂直的網(wǎng)格切片被給出;計算域(邊界條件已在圖中標出)的尺寸為10C(x)×4C(y)×6C(z)。在胸鰭附近的核心區(qū)域內(nèi)以間距為0.008 6C的均勻網(wǎng)格提供極高的分辨率。在核心區(qū)域之外,網(wǎng)格向計算域的外邊界迅速拉伸(對于近尾流區(qū)保證拉伸因子在3%以內(nèi))。網(wǎng)格無關(guān)性研究通過將胸鰭表面網(wǎng)格數(shù)和核心區(qū)域內(nèi)3個坐標方向上的網(wǎng)格數(shù)加倍來實施,同時為保證整個網(wǎng)格系統(tǒng)的光滑性外部區(qū)域的網(wǎng)格數(shù)也做適當增加,細化后的流場網(wǎng)格總數(shù)為2 355萬。計算域無關(guān)性研究是將計算域尺寸在3個坐標方向上擴大一倍來進行,在此過程中核心區(qū)域及其網(wǎng)格劃分保持不變,而在外部區(qū)域使用與上述標準網(wǎng)格同樣的拉伸因子。
圖4給出了在數(shù)值計算中最后3個周期內(nèi)(4T~7T)推力系數(shù)隨時間的變化情況,此時流動已達周期性穩(wěn)定狀態(tài)(時間步長為T/400)??梢园l(fā)現(xiàn),在標準網(wǎng)格(nominal grid)、細網(wǎng)格(finer grid)和更大的計算域(larger domain)上計算得到的推力系數(shù)非常接近。因此,在接下來的模擬中采用標準網(wǎng)格,這樣既能保證計算結(jié)果的精確性,也能使計算所需的時間和內(nèi)存減少。
圖5給出了在相同條件下的試驗結(jié)果和各種數(shù)值方法的計算結(jié)果。經(jīng)過7個運動周期的模擬流場已達周期性穩(wěn)定狀態(tài),圖中的是最后一個運動周期的模擬結(jié)果??傮w來看,浸入邊界法和Wang等[23]的試驗結(jié)果,而蘇玉民等[15]的勢流理論方法由于沒有考慮流體的黏性而與試驗結(jié)果差距較大。進一步地對比可以發(fā)現(xiàn),浸入邊界法計算得到的CT隨無量綱時間的變化曲線有與試驗吻合較好的平緩波峰,而從FLUENT計算結(jié)果中則可看到較為尖銳的波峰。
圖3 計算域和網(wǎng)格分布Fig. 3 Computational domain and grid distribution
圖4 不同網(wǎng)格和計算域下推力系數(shù)隨時間的變化Fig. 4 Variation of thrust coefficient with time for different grids and computational domains
圖5 推力系數(shù)數(shù)值與試驗結(jié)果比較Fig. 5 Comparison of numerical and experimental results of thrust coefficient
3.2.1 推力系數(shù)和升力系數(shù)的脈動
圖6給出了一個運動周期內(nèi)胸鰭的水動力系數(shù)隨時間的變化。其中雷諾數(shù)和斯特勞哈爾數(shù)的取值分別為Re=2 500和St=0.55,所選取的Re和St值均包含在Lauder和Jayne[17]的黑鱸魚胸鰭試驗所給出的無因次參數(shù)范圍(1 756≤Re≤4 389和0.32≤St≤0.89)之內(nèi),還將在3.3節(jié)中對St的影響進行討論。數(shù)值計算中其他的運動參數(shù)為φRC=30°,φFEC=-30°,φRA=30°,φFEA=30°,ΔφEF=90°,類似的參數(shù)取值已在有關(guān)胸鰭推進性能分析的有關(guān)文獻[9, 18]中被采用。從圖6中可以看到,胸鰭推力系數(shù)和升力系數(shù)曲線在每個運動周期均只有一個峰值。根據(jù)當前的計算,胸鰭在前半個周期主要產(chǎn)生推力,而在后半個周期主要產(chǎn)生升力。鑒于此在有關(guān)試驗和數(shù)值研究[24]中,通常把胸鰭在一個周期內(nèi)的劃水過程分成兩個階段:動力劃水(power stroke)階段和恢復(fù)劃水(recovery stroke)階段。
圖6 水動力系數(shù)隨時間的變化Fig. 6 Variation of hydrodynamic coefficients with time
3.2.2 尾渦結(jié)構(gòu)與推進機理分析
一個周期內(nèi)胸鰭運動的8個典型瞬時的三維尾渦結(jié)構(gòu)(左圖)、下游表面壓力分布(中間圖)和上游表面(迎流面)壓力分布(右圖)如圖7所示。三維空間渦結(jié)構(gòu)的識別采用Q判據(jù)[25],左圖中渦結(jié)構(gòu)的顏色反映了渦量的大小,兩側(cè)表面壓力分布圖中胸鰭上的矢量為局部表面力。
1) 動力劃水階段
圖7(a)給出了一個運動周期的初始時刻的情況。梢渦TV1、腹側(cè)邊緣渦VEV1和背側(cè)邊緣渦DEV1已從仿生鰭上脫落,尾流場中的兩個渦環(huán)R1和R2為前面的兩個運動周期所產(chǎn)生。
再來看1/8周期時的情況,如圖7(b)所示,胸鰭上游表面有大片的低壓力區(qū)出現(xiàn),這是由附著在鰭邊緣的馬蹄形渦結(jié)構(gòu)所誘導的,該渦結(jié)構(gòu)的形成與胸鰭的快速向后擺動有關(guān),它由新梢渦TV、新腹側(cè)邊緣渦VEV和新背側(cè)邊緣渦DEV連接而成。大片的低壓力區(qū)意味著該運動瞬時仿生鰭的上游和下游表面之間有大壓差力存在,而且該瞬時胸鰭表面和上游來流之間好的垂直度也有利于壓差力在推力方向的分配,所以仿生鰭的推力在該瞬時幾乎達到峰值(參考圖6(a))。
圖7 典型瞬時的三維尾渦結(jié)構(gòu)和鰭表面壓力分布Fig. 7 Three-dimensional wake structure and pressure distribution on the fin surface at typical instants
圖7(c)給出了1/4運動周期時的情況,可以看出,梢渦TV和背側(cè)邊緣渦DEV從鰭上脫落的過程已經(jīng)開始,R1和R2兩個渦環(huán)伴隨著一定程度的耗散被尾流帶向下游。
下一個運動相位(t/T=3/8)如圖7(d)所示。TV連同部分的DEV和部分的VEV從胸鰭上脫落,這導致迎流面的低壓區(qū)相比于前一時刻(t/T=1/4)繼續(xù)減小。
2) 恢復(fù)劃水階段
在圖7(e)中,即t/T=1/2的運動相位下,渦結(jié)構(gòu)TV+VEV+DEV已經(jīng)完全從胸鰭上脫落,鰭梢處又有新梢渦TVW正逐漸成長。
圖8 胸鰭特有的三維雙環(huán)渦Fig. 8 Particular three-dimensional dual-ring vortex for pectoral fin
t/T=5/8運動瞬時的情況如圖7(f)所示,圖7(e)中脫落的馬蹄形渦與TVW相連而產(chǎn)生了一個新的渦環(huán),胸鰭的腹側(cè)邊緣和背側(cè)邊緣又分別有新的附著渦VEVW和DEVW生成。
圖7(g)給出了t/T=3/4時的尾渦結(jié)構(gòu)和壓力分布,新梢渦正在RW的誘導下被拉伸。仿生鰭向前擺動的同時,DEVW和VEVW上的渦量不斷積累。在強大的DEVW的誘導下,一個顯著的低壓力區(qū)在胸鰭的背側(cè)邊緣附近出現(xiàn),同時下游表面該區(qū)域存在較大壓差力的事實也可從表面力矢量分布中看出,而且在這個運動相位時鰭面與升力軸之間的垂直度較好,導致被分解到升力方向的壓差力比重增加。在上述這些因素的共同作用下,仿生鰭在3/4運動周期時的升力接近峰值(參考圖6(b))。
圖7(h)給出了7/8周期時的情況,此時仿生鰭在本周期的劃水過程即將完成。TVW、VEVW和DEVW正在從鰭上逐步脫落,下一個運動周期的TV1、VEV1和DEV1(見圖7(a))實際上正是從這些渦結(jié)構(gòu)發(fā)展而來,同時RW這個渦環(huán)將發(fā)展成圖7(a)尾流場中的渦環(huán)R1,這樣下一個周期的仿生鰭運動和三維尾渦演變的循環(huán)即將開始。
3) 胸鰭尾渦結(jié)構(gòu)的關(guān)鍵特征——三維雙環(huán)渦結(jié)構(gòu)。
以上關(guān)于推進機理的探討事實上揭示出胸鰭三維尾渦的一個重要特征——雙環(huán)渦結(jié)構(gòu),該特殊渦結(jié)構(gòu)是由動力劃水所生成渦環(huán)RW(R1)與恢復(fù)劃水所生成渦環(huán)TVW+VEVW+DEVW (TV1+VEV1+DEV1)以接近垂直的角度相連而形成的。為使這一點變得更加清晰,繪制了圖8,其中Rpower表示動力劃水階段所產(chǎn)生的渦環(huán),它與恢復(fù)劃水渦環(huán)Rrecovery一起構(gòu)成了胸鰭獨特的雙環(huán)渦結(jié)構(gòu)DR。
對于仿生翼/鰭推進及魚類游動來說,斯特勞哈爾數(shù)St是一個關(guān)鍵的無因次數(shù),它特征化了流場渦結(jié)構(gòu)和水動力性能[7]。就仿生翼及魚類尾鰭而言,研究表明推進效率在St=0.3附近取得最大值[7, 26-27]。本小節(jié)分析St對胸鰭尾流場和推進性能的影響,及考查胸鰭推進的最優(yōu)St。
1) 三維尾渦結(jié)構(gòu)
在胸鰭一個劃水動作的循環(huán)完成(t/T=1.0)時,不同St下的三維流場渦結(jié)構(gòu)如圖9所示,其他運動參數(shù)的取值與3.2節(jié)一致。能夠觀察到,胸鰭運動產(chǎn)生的尾渦隨斯特勞哈爾數(shù)的增大而在渦強上有所增加,特別是對于胸鰭附近的渦結(jié)構(gòu)(如TVW+VEVW+DEVW)而言。另外,模擬也發(fā)現(xiàn)在低St(0.25)時,前半個周期所生成的渦環(huán)RW與后半個周期所生成的渦環(huán)TVW+VEVW+DEVW分離。同時也可以看到,遠尾流區(qū)的渦環(huán)取向隨斯特勞哈爾數(shù)的增加而有所改變(即水平軸與渦環(huán)軸向間的夾角變小),從而在渦環(huán)所引起的射流中流向動量成分的比例提升。
圖9 尾渦結(jié)構(gòu)隨斯特勞哈爾數(shù)的變化Fig. 9 Variation of wake structure with Strouhal number
2) 水動力性能
St的改變對推力系數(shù)和效率的影響如圖10所示。可以看出,胸鰭的平均推力隨斯特勞哈爾數(shù)的增大而增加,這是由流場渦結(jié)構(gòu)強度的增加和渦環(huán)的取向變化所導致的(見圖9),這兩個因素造成與推力產(chǎn)生密切相關(guān)的射流的流向動量變大,進而引起推力的增加。
另一方面,效率隨St的增加而先增大后減小。定量來看,效率先從22%(St=0.25)較快地提升到44%(St=0.55),然后又較緩地降低到32%(St=0.85)。開始時效率隨St的提升可解釋為推力隨斯特勞哈爾數(shù)的增大。達到最佳后效率隨St的降低有以下兩個原因: 1) 尾渦的強度增加(參考圖9)使得更多的能量被尾流所攜帶; 2) 各渦結(jié)構(gòu)之間的聯(lián)系增多(參考圖9)使得黏性抵消效應(yīng)變得更顯著,從而引起胸鰭向后射流動能的額外損耗。
以上數(shù)值模擬及分析表明:就胸鰭推進來說,也同樣存在最優(yōu)的St使效率最佳。本節(jié)的研究指出胸鰭的最優(yōu)斯特勞哈爾數(shù)約為0.55,這與Lauder等[11]的試驗得到的胸鰭進行常規(guī)游動時的斯特勞哈爾數(shù)的值0.54很接近。
圖10 不同斯特勞哈爾數(shù)下的水動力性能Fig. 10 Hydrodynamic performance for different Strouhal numbers
1) 三維尾渦結(jié)構(gòu)
在胸鰭一個劃水動作的循環(huán)完成(t/T=1.0)時,不同St下的三維流場渦結(jié)構(gòu)如圖11所示,其他運動參數(shù)的取值與3.2節(jié)一致。不難發(fā)現(xiàn),在所研究的相位差范圍內(nèi)胸鰭的尾渦結(jié)構(gòu)是相似的,但ΔφFE=45°時的鰭梢部渦結(jié)構(gòu)TVW相比于兩個更高的相位差的情況更強。也可以看出,ΔφFE=45°時胸鰭的尾流寬而短,這種發(fā)散的尾流通常是低推力和效率的指示。在高相位差135°時,一個值得注意的現(xiàn)象是特殊柱狀渦HV的出現(xiàn),它由呈螺旋形的渦絲構(gòu)成,在升力翼的尾流場中也出現(xiàn)了類似的螺旋形渦結(jié)構(gòu)[28],比較90°和135°時的平均升力模擬結(jié)果,發(fā)現(xiàn)后者比前者高了21%。Bozkurttas等[12]在胸鰭尾流場中也觀察到了相似的螺旋形柱狀渦。
圖11 尾渦結(jié)構(gòu)隨相位差的變化Fig. 11 Variation of wake structure with phase angle
2) 水動力性能
圖12給出了推力系數(shù)和效率隨相位差的變化。當前的計算表明,在ΔφFE=90°時胸鰭同時取得最高的推進效率(44%)和平均推力系數(shù)(4.91)。當相位差以90°為標準值增加或減小時,均可以看到推力和效率的降低。具體地,相位差為135°和45°時的平均推力系數(shù)分別為3.57和2.23,推進效率分別為24%和15%。胸鰭在低相位差(45°)時產(chǎn)生較小的推力和效率與該相位差下的發(fā)散尾流(見圖11)有關(guān),而且從定量計算結(jié)果來看相位差以90°為中心上升和下降時所引起的胸鰭推進性能指標下降并不是對稱的,胸鰭在高相位差(135°)時效率的降低與升力增加(相比90°時增加了21%)所導致的在垂向方向的能量損失增多有關(guān)。
圖12 不同相位差下的水動力性能Fig. 12 Hydrodynamic performance for different phase angles
通過仿生胸鰭的三維尾渦結(jié)構(gòu)與參數(shù)影響分析,得到如下結(jié)論:
1) 在一個周期中,胸鰭的劃水動作由主要產(chǎn)生推力的動力劃水階段和主要產(chǎn)生升力的恢復(fù)劃水階段構(gòu)成。推進機理研究揭示出胸鰭的尾流場由一個特殊的三維雙環(huán)渦所支配,其中動力劃水渦環(huán)誘導更多的流向動量成分,而恢復(fù)劃水渦環(huán)誘導更多的垂向動量成分。
2) 隨著St的增大,胸鰭的推力單調(diào)增加,而效率經(jīng)歷了一個先較快提升后較慢下降的過程。背后的機理是尾渦結(jié)構(gòu)的強度隨St的增加而增大,渦環(huán)中心軸的方向與來流方向之間的夾角減小,同時尾流所帶走的能量增多,各渦結(jié)構(gòu)之間的相互作用增強。
3) 胸鰭推進的最優(yōu)St范圍比拍動翼和尾鰭推進的更高。當配有仿生胸鰭的機器魚進行常規(guī)游動時更高的效率是我們所希望的,因此胸鰭的St應(yīng)處于上述最優(yōu)范圍;而當機器魚加速起動時更大的推力是我們所希望的,因此胸鰭應(yīng)該工作在更高的St下以產(chǎn)生高推力,同時模擬也表明高St下胸鰭仍具有可觀的效率。
4) 胸鰭的推力和效率達到最大值時的相位差均為90°。低相位差時寬而短的發(fā)散尾流造成了較低的推進性能,而高相位差時的螺旋形柱狀渦導致了較高的升力。
5) 當配有仿生胸鰭的機器魚進行常規(guī)游動時,相位差應(yīng)保持在90°以取得最佳的推力和效率;而當機器魚進行操縱運動時,胸鰭應(yīng)工作在更大的相位差下以獲得更多的升力。