黃愷俊 余永亮

(中國科學院大學生物運動力學實驗室,北京 100049)

引言

自然界中的大多數(shù)魚類進化出波狀擺動的方式實現(xiàn)推進[1-3],其中振幅包絡(luò)線從前部到后部逐漸增加的波狀擺動更有利于產(chǎn)生推力[4],這類波動模式通常被稱為鲹科模式.鲹科模式與鰻鱺模式、金槍魚模式都是通過身體和尾鰭的擺動進行推進的,鲹科模式的魚體前 1/3不波動,后體擺動的波幅隨體長增大[3],常見的斑馬魚、草魚、鯉魚等都采用這類推進模式.鰻鱺模式的魚類全部魚體都大幅擺動,而金槍魚模式的魚體僅尾部往復(fù)擺動.鲹科魚類的游速較高,常被用于水下仿生航行器的仿生設(shè)計中.因此,深入研究鲹科模式推進的流體力學機制不僅可以增強對波狀推進的認識,還可以為水下航行器的設(shè)計提供參考.

Videler[2]和Lighthill[5]預(yù)測了魚尾擺動后的流場尾跡,描述波狀推進的魚類在尾部會脫落出交錯排列的離散尾渦,被稱為反卡門渦街.前人通過數(shù)字粒子圖像測速儀(DPIV)觀測到太陽魚的尾跡中存在交替的渦結(jié)構(gòu)[6-8],在振蕩翼的模型實驗中也發(fā)現(xiàn)相同的流動[9-11];通過對二維波狀翼型的數(shù)值模擬,Dong 等[12]、Deng 等[13]和Liu 等[14]也同樣觀察到這種典型的尾跡結(jié)構(gòu).與阻力型的卡門渦街不同,反卡門渦街是一種推力型的尾流結(jié)構(gòu).Yu 等[15]基于波狀推進的魚體產(chǎn)生反卡門渦街的共性,利用渦動力學分析方法,推導并提出了統(tǒng)一描述波狀推進的水動力與運動參數(shù)、雷諾數(shù)的關(guān)聯(lián)關(guān)系的標度律.該規(guī)律對鰻鱺模式、鲹科模式以及金槍魚模式都適用.

在研究魚體運動與受力關(guān)系方面,近10 年來也有不少實驗和計算分析的工作.Gazzola 等[16]結(jié)合實驗和數(shù)值模擬的結(jié)果,通過流體力中推、阻力的匹配關(guān)系,得到巡游速度的標度律.Thekkethil 等[17]使用統(tǒng)一的運動學模型進行二維數(shù)值模擬,該模型涵蓋各種魚體波動和魚尾剛體俯仰運動,并研究了波動和俯仰運動的等價性以及推力產(chǎn)生的機制.Hu 等[18]對魚體的低頻波狀擺動進行數(shù)值模擬,獲得魚體阻力與波動參數(shù)之間的標度律.Gupta 等[19]則不光給出波動魚體受力,還研究了輸出功率、效率的標度律關(guān)系.Gao 等[20]在Yu 等[15]的工作基礎(chǔ)上,基于壓差力標度律和摩阻力標度律的匹配關(guān)系,也給出巡游速度的標度律.這些工作很好地給出運動與水動力性能的宏觀聯(lián)系,為進一步分析推進的流體力學機制提供數(shù)據(jù)支持.

事實上,通過計算手段獲得流場數(shù)據(jù)只是全面了解流動的第一步,為了找到流動和受力的定量關(guān)系,引入氣動力理論的方法可以從流動數(shù)據(jù)中提取物理規(guī)律[21].目前有兩種方法可以建立起物體所受力與流場結(jié)構(gòu)之間的聯(lián)系.一種是Burgers[22],Lighthill[23]和Wu[24]提出的導數(shù)矩的方法,通過流體的渦量一階矩計算流場動量,并通過動量變化求得物體受力.后來被Wu 等[25]發(fā)展到有限域的渦量矩理論.另一種是Quartapelle 等[26]、Howe[27]和Chang[28]采用的投影法,通過引入一個勢函數(shù),將力分解為附加質(zhì)量力,蘭姆矢量相關(guān)的旋渦貢獻,以及黏性相關(guān)的項.Yu[29]引入虛速度的概念,通過變分得到了流體力學中的虛功率原理,并指出引入的勢函數(shù)實際上是物體虛運動引起的虛速度場,虛速度場具有遠場衰減的特性,在受力分解中包含物面積分項和流體結(jié)構(gòu)的體積分項,并分析各項的伽利略不變形式.為了研究鲹科模式波狀推進的流體動力響應(yīng),從虛功率原理出發(fā),可以深入分析變形體的邊界瞬時流體力效應(yīng)和流場結(jié)構(gòu)對流體力的貢獻.

本文選取鲹科模式推進的魚體作為研究對象,研究魚游波狀推進的流體力機理.第1 節(jié)介紹魚體波動模型以及基于虛功率原理的分析方法.第2 節(jié)通過研究典型參數(shù)下鲹科魚波狀擺動流體力的分解來闡明推力的主要來源,然后討論不同雷諾數(shù)條件下推力與Strouhal 數(shù)的關(guān)系,并結(jié)合波狀推進的標度律分析,研究標度律公式中各項物理意義.

1 計算模型與分析方法

1.1 魚游模型

借鑒前人[12,17,19,30]研究鲹科模式魚類游泳的物理模型,本文采用NACA0012 翼型來模擬魚體,翼型的中線代表了魚體在靜態(tài)平衡時的脊柱.魚體脊柱波動的側(cè)向位移y(x,t)由Videler 等[31]通過實驗提供的行波方程來模擬

其中,λ和c分別表示波動波長和波速,波速是頻率f與波長 λ的乘積.如圖1 所示,魚體靜止時,頭部位于x=0,尾部在x=L處.Am(x)表示不同身體位置的最大振幅.

圖1 仿魚翼型波狀擺動示意圖Fig.1 Schematic of a fish-like undulatory airfoil

對于鲹科模式的波動,通常認為波長 λ等于一倍體長,身體的振幅包絡(luò)線Am(x)可由二次函數(shù)給出

其中Am(0)=0.02L,Am(0.2)=0.01L,Am(1)=0.1L.圖2給出包絡(luò)線示意圖(虛線)和魚體中線的側(cè)向位移波動曲線(實線),可以看出,最大振幅出現(xiàn)在尾尖處,而最小振幅大概在 0.23L處.魚體頭部有小幅的波動,而尾部相對于頭部有更高的振幅.本文定義魚體尾部最大振幅為魚體波動的振幅,即波動振幅A=0.1L.

圖2 一個波動周期不同時刻的魚體中心線及其包絡(luò)線Fig.2 The centerline of the fish body at various time instants for one undulation cycle and their envelops

在本文中,由波動魚模型引起的二維非定常流動,可由以下無量綱化的不可壓縮Navier-Stokes (N-S)方程控制

其中,u是速度場,p是壓強場,μ為流體介質(zhì)的黏性系數(shù).本文選取特征的參考速度U(來流速度),魚體的體長L以及流體密度 ρ對所有的物理量進行無量綱化,則體現(xiàn)流場慣性力和黏性力比值大小的雷諾數(shù)為Re=ρUL/μ.另一個常用的無量綱數(shù)是Strouhal數(shù),它表征了流動的非定常性,定義為St=2f A/U.本文采用的數(shù)值求解N-S 方程的程序是開源程序包OpenFOAM 的pimpleFoam 求解器,程序的驗證和確認見附錄,該程序在文獻[15]中已得到了驗證.

1.2 基于虛功率原理的推力分析方法

通過數(shù)值手段獲得魚體波狀擺動的流場信息,便可通過魚體表面的法向應(yīng)力pn和切向應(yīng)力 τ=μn×ω積分,得到魚體的受力F

其中,法向量n指流體的外法向量,和魚體表面的外法向量方向相反.該公式簡潔地給出魚體與流體的相互作用結(jié)果,可用于驗證以下基于虛功率原理的推力分析方法.

為了建立流動與受力的定量聯(lián)系,本文引入Yu[29]提出的虛功率原理的分析方法.假定魚體在前進方向上產(chǎn)生瞬時的虛運動,則流場中的虛功率關(guān)系可得

其中,u*是由魚體虛運動引起的虛速度場,σji代表施加在流體微團上的表面力的張量(稱為應(yīng)力張量).式(5)的物理意義是物體所受水動力沿著虛運動方向所做的虛功率等于流場中的內(nèi)應(yīng)力和廣義力的虛功率.廣義力包含流場中的體積力fi和慣性力 dui/dt.進一步地,方程(5)左側(cè)的值就是推力FT,即水動力沿著虛運動方向的分量,所以對方程(5)右側(cè)各項的分析將有助于理解推力的物理成因.

Yu[29]在文中提到,可以假定虛流場是不可壓縮勢流.為了使勢函數(shù)唯一,采用了在無窮遠處勢函數(shù)φ為0 的邊界條件[26,28-29].因此,虛流場滿足的方程和邊界條件如下

其中 ? φ是虛流場的速度矢量u*,Sfish和S∞分別表示魚體表面和無窮遠邊界.虛流場是無源無旋的,可通過在魚體表面分布基本解(源匯和偶極子)以滿足方程(6)的魚體和無窮遠處的邊界條件.將虛流場方程(6)代入虛功率關(guān)系(5),并考慮流動不可壓縮且體積力可忽略,則表達式(5)右式的兩項分別為

從式(7)可以看出流場的流體內(nèi)應(yīng)力對虛應(yīng)變做功的功率轉(zhuǎn)換為壁面的切應(yīng)力 τ所做的虛功率,一部分是對虛速度場做功的功率,另一部分是對魚體虛運動做功的功率.前者被稱為類摩阻(該說法來自Chang[28]),后者正好為摩阻在推進方向的大小.從式(8)可以看出流場中慣性力的虛功率可以轉(zhuǎn)換為與魚體邊界加速度有關(guān)的表面積分和與流場中對流散度有關(guān)的體積分.表面積分表明壁面的瞬時運動即可確定部分虛功率的大小,其值可以看作在黏性流動中魚體運動狀態(tài)瞬時所確定的推力貢獻.而流體對流的散度正好等于負的2 倍Q(Q是速度梯度張量的第二不變量,定義為Q=0.5(||Ω||2-||D||2),其中 Ω 和D分別是速度梯度分解成的反對稱張量和對稱張量),Q＞0代表流體質(zhì)點的旋轉(zhuǎn)強于應(yīng)變速率,Hunt 等[32]就提出用Q＞0作為旋渦的識別.因此,式(8)中右邊第二項可以看作是流場中流體的旋轉(zhuǎn)和應(yīng)變速率相對大小對推力的貢獻,為了方便,我們稱之為Q貢獻.

因此,本文將魚體受到的推力FT分解為4 個分量

分別表示邊界加速度對推力的貢獻FT,a、流場中Q對推力的貢獻FT,Q、邊界的類摩阻分量FT,v以及摩阻分量FT,f,其推力系數(shù)及各分量的系數(shù)定義為

2 結(jié)果與討論

2.1 典型問題的推力分解

為了考察式(9)受力分解與式(4)在推力分析上的一致性,選擇無量綱擺尾幅度為A=0.10,波長為 λ=1.0,波速為c=2.0,流動的雷諾數(shù)為Re=5000

作為典型的案例,對比了一個擺動周期里推力系數(shù)隨時間的變化,如圖3(a)所示.結(jié)果表明基于虛功率原理的推力曲線與由式(4)直接積分的推力曲線重合,具有一致性.進一步地,圖3(b)給出了4 個分量CT,a,CT,Q,CT,v和CT,f隨著時間變化的曲線,發(fā)現(xiàn)邊界加速度對推力始終是一個正的貢獻,而流場中Q對推力有很大的負貢獻,兩者相位幾乎相反.邊界的類摩阻分量相對于其他分量要小得多,幾乎可以忽略不計,而摩阻分量在一個周期內(nèi)始終是阻力.四者的周期均值分別為 0.2474,-0.1192,0.0003和-0.0607,可見CT,v的確可以忽略,推力主要由CT,a,CT,Q,CT,f3 項組成,其中CT,Q接近摩阻CT,f的兩倍.接下來我們將對推力貢獻最重要的兩部分CT,a和CT,Q做進一步的分析.

圖3 推力系數(shù)在一個周期內(nèi)的變化Fig.3 The time-dependent thrust coefficient in an undulation cycle

對于二維波動模型,可以對邊界加速度對推力的貢獻FT,a做進一步的推導.分別沿著上、下表面進行積分可以得到

其中,下標+和下標-分別代表上下表面的物理量,y+和y-分別代表上下表面的曲線方程,y′=?y/?x代表曲線對x的導數(shù).根據(jù)變形規(guī)則,在相同的x坐標,魚體邊界與魚體脊柱波動一致,則邊界加速度的值du/dt是相同的.邊界加速度可由波動方程對時間的二階導數(shù)獲得,流體在魚體表面的外法向量可由曲線方程獲得,具體關(guān)系如下

將方程(12)代入方程(11),我們可以得到

基于式(13)我們可以看出由邊界加速度引起的推力貢獻沿著魚體中線的分布與兩部分相關(guān)因素相關(guān).一部分是側(cè)向運動的y方向加速度(波動加速度),另一部分是上下表面的虛速度勢的差值.根據(jù)勢流的特點,魚體表面的虛速度勢是由當前時刻的魚體外形決定的,因此

圖4(a)給出了一個擺尾周期內(nèi)邊界加速度對推力貢獻沿著中線的分布,其中實線和點被用來區(qū)分擺尾前半周期和后半周期的特定時刻.從圖4 中可以看出,前半周期和后半周期的擺尾動作產(chǎn)生的FT,a是一致的.從曲線看,頭部區(qū)域由于變形小,加速度小,對FT,a的影響就小,尾部區(qū)域的值波動大,大部分時刻都處于正值,這意味著對推力有大貢獻.進一步地,為了更好地展示魚體各個部位變速運動帶來的貢獻,將魚體沿中線等分成5 個區(qū)域,并統(tǒng)計了5 個區(qū)域分別對推力貢獻的周期均值,如圖4(b)所示.從數(shù)據(jù)上可以看出,前3 個區(qū)域的均值與整體均值之比大約為3%～4%,且均為負貢獻,處于尾部的第4 和第5 區(qū)對整體正推力有很大的正貢獻,特別是加速度最大的第5 區(qū),雖然只占魚體長度的20%,但這個區(qū)域的貢獻占整個貢獻的83.25%.

圖4 邊界加速度對推力的貢獻Fig.4 The contribution of the boundary acceleration

根據(jù)式(9)右式第二項可知,流場中的Q分布對推力的貢獻依賴兩個因素,一個是Q的強度,另一個是虛速度勢 φ.為了更好地展現(xiàn)流場信息,圖5 給出了4 個典型時刻的流場渦量分布、Q值分布、Q對推力貢獻分布,以及虛速度勢的分布.從渦量場來看,魚體上半部分邊界層內(nèi)為負渦量,下半部分為正渦量,尾跡呈現(xiàn)典型的反卡門渦街結(jié)構(gòu).從Q值圖上可以看出,在魚體上下部分的邊界層處,Q值僅在尾部附近有明顯的分布值,且兩側(cè)的值相反,意味著尾部拍動的過程中,“迎風面”(圖中顯示的魚體尾部下表面)邊界層內(nèi)的流體應(yīng)變速率強于旋轉(zhuǎn)(Q＜0),而另一側(cè)的尾部附近流體的旋轉(zhuǎn)強于應(yīng)變速率(Q＞0).在魚體頭部,存在很強的應(yīng)變速率區(qū)(Q＜0),頭部兩側(cè)后方存在很強的旋轉(zhuǎn)區(qū)(Q＞0).對于脫落到流場中的自由旋渦,無論是正渦還是負渦,中心都由流體旋轉(zhuǎn)主導(Q＞0),且在旋渦周圍存在較弱的應(yīng)變速率區(qū)(Q＜0).隨著魚體的擺動,虛流場的速度勢函數(shù) φ也會隨之變化,特別是在尾部區(qū)域,魚體外形變化大導致 φ值變化大.流場中流體的旋轉(zhuǎn)和應(yīng)變速率的相對大小對推力的貢獻(-2ρQφ) 在頭部區(qū)域有正有負,在尾部邊界層區(qū)域也有正有負.在尾跡區(qū),由于虛速度勢始終為正(φ ＞0),故脫落出去的旋渦中心貢獻了負推力,旋渦周圍呈現(xiàn)一個較弱的正貢獻.

圖5 t=0T,t=T/8,t=T/4,t=3T/8 時刻的渦量場(第1 列)、Q 分布圖(第2 列)、Q 對推力的貢獻分布圖(第3 列)和虛速度勢場(第4 列)Fig.5 Vorticity field (first column),Q distribution diagram (second column),Q distribution diagram of contribution to thrust (third column) andvirtual velocity potential field φ (fourth column) at t=0T,t=T/8,t=T/4,t=3T/8

為了更好地定量展示流動結(jié)構(gòu)對推力的貢獻,我們基于流場中Q分布的特點,人為對整個流場進行了分區(qū).如圖6 所示,zone0 代表了距離魚頭部一定距離的流場,該區(qū)域的Q值都很小.zone1 是魚體頭部附近的區(qū)域,其Q分布有正有負.zone2～zone5 分別代表身體波動的部分,且波動的振幅逐漸增強,這些區(qū)域分布的Q強度也逐漸增強.區(qū)域zone_wake 則是尾跡區(qū),該區(qū)域又被分為對應(yīng)的5 個子區(qū),由近及遠分別為wake1～wake5.

圖6 流場分區(qū)示意圖Fig.6 Schematic of flow field zoning

圖7(a)給出了各個分區(qū)Q分布對推力的貢獻隨著時間變化的曲線,圖7(b)給出了各個分區(qū)對推力貢獻的周期均值,它們具有如下幾個特點.(1)從時變曲線來看,zone0 和zone1 的推力貢獻值都在零附近小幅振蕩.由于zone0 區(qū)的流場擾動小,Q值小,整體流動幾乎沒有推力貢獻.zone1 的頭部區(qū)存在強正強負的貢獻,但整個擺動周期中變化小,其貢獻也較小.(2)魚體波動的zone2～zone5 這4 個區(qū)因擺動振幅逐漸增大而對推力貢獻值也波動增強.其中,zone2,zone3 和zone4 的波動曲線仍在零值附近振蕩,從均值來看,zone2 和zone3 區(qū)域?qū)ν屏τ胸暙I,但zone4 區(qū)呈現(xiàn)出阻力特性.zone5 是魚體尾尖附近的區(qū)域,振蕩最強,其周期均值最大,具有大阻力特征.(3)尾跡區(qū)在擺動中大部分時間里都呈較弱的振蕩,其周期均值為 9.31%,體現(xiàn)的是阻力,遠小于zone5的作用.由此可見,流場中的流動結(jié)構(gòu)帶來的阻力振蕩也主要來自于尾部附近.

圖7 Q 分布對推力的貢獻Fig.7 The contribution of Q distribution to thrust

為了更好地展示尾跡區(qū)的旋渦貢獻,圖8 給出了尾流區(qū)(wake1～wake5)的Q對推力貢獻隨時間的變化.結(jié)果表明,每一個分區(qū)對推力的貢獻都是負的,越遠離魚尾的尾跡區(qū)域,曲線的漲落就越小,均值也越小.曲線的漲落主要是由于旋渦在區(qū)域內(nèi)流進流出引起的,而遠離魚體的尾跡區(qū)的均值減小是由于虛速度勢函數(shù)是隨距離的衰減函數(shù).整體來看,尾跡區(qū)對魚體受力的影響是較弱的.

圖8 不同尾跡區(qū)的 Q 分布對推力的貢獻Fig.8 The time-dependent contribution of Q distribution to thrust in different wake regions

本小節(jié)討論了典型波狀推進情況下的流體力分解.研究發(fā)現(xiàn),魚體波狀擺動產(chǎn)生的推進力,其正貢獻主要來自于邊界加速度的瞬時響應(yīng),而流場中的流體旋轉(zhuǎn)和應(yīng)變速率相對大小對應(yīng)的是阻力,在鲹科類的擺尾推進中,兩者都主要體現(xiàn)在尾部附近的流動.除此之外,邊界上的黏性耗散效應(yīng)主要體現(xiàn)在摩阻上.

2.2 推力與St 數(shù)關(guān)系

圖9 不同雷諾數(shù)Re條件下,〈CT,a〉,〈CT,Q〉,〈CT,v〉和〈CT,f〉 隨運動參數(shù)St的變化規(guī)律Fig.9 The time-averaged thrust coefficient 〈CT,a〉,〈CT,Q〉,〈CT,v〉,〈CT,f〉 as a function of St at different Reynolds number

為了更好地顯示4 部分時均值對推力貢獻的相對大小,圖10 給出了Re=500 和Re=5000情況下的各項貢獻的堆疊柱狀圖.從圖中可以看出,無論是低雷諾數(shù)還是高雷諾數(shù)情況,正推力的貢獻都主要來自于邊界變速運動.壁面切應(yīng)力的類摩阻分量無論產(chǎn)生正貢獻還是負貢獻,相對于其他3 項總是較小的.阻力分量主要來源于流場中Q的貢獻和壁面摩阻的貢獻.對于低雷諾數(shù)(Re=500)的情況,兩者的阻力貢獻相當,摩阻的阻力貢獻略大于流場中Q的貢獻.對于高雷諾數(shù)(Re=5000) 的情況,流體中Q分布的阻力貢獻強于摩阻的貢獻.

圖10 不同雷諾數(shù)Re條件下,〈CT,a〉,〈CT,Q〉,〈CT,v〉和〈CT,f〉的貢獻堆疊柱狀圖Fig.10 Stacked histogram of the time-averaged thrust coefficient 〈CT,a〉,〈CT,Q〉,〈CT,v〉,〈CT,f〉 at different Reynolds number

圖11 在 t=0T,t=T/8,t=T/4和t=3T/8 時刻以及不同擺幅條件下(A=0.06,0.08,0.10,0.12,0.14),Δφ(x) 以及 Δφ(x)/A 沿著魚體中線的分布Fig.11 Distribution of Δφ(x) and Δφ(x)/A along the center line of fish at different time (t=0T,t=T/8,t=T/4,t=3T/8) and under different amplitude condition (A=0.06,0.08,0.10,0.12,0.14)

曲線幾乎重合,這表明 Δφ(x)與擺幅A成正比,即,Δφ(x)∝A.

因此,對于任意時刻,FT,a沿著魚體中線的分布都與A2f2成正比,由此可得

即由邊界加速度得到的瞬時推力和推力均值都與S t2成正比,并且圖4(b)給出的分區(qū)貢獻占比不會隨著運動參數(shù)S t的變化而變化.

進一步關(guān)于時均值 〈CT,a〉,〈CT,Q〉,〈CT,v〉和〈CT,f〉的分析可以結(jié)合已經(jīng)發(fā)表的關(guān)于魚體波狀推進標度律[15]進行分析.

2.3 波狀推進的標度律分析

Yu 等[15]基于渦動力學分析方法,給出了仿魚類波狀擺動推進器的推力系數(shù)與S t和Re之間的標度關(guān)系,即

圖12 慣性力對推力貢獻的標度關(guān)系及其擬合參數(shù)與雷諾數(shù)的關(guān)系Fig.12 Scaling law of the contribution of the inertial force to thrust and the relationship between fitting parameters and Reynolds number

圖13 內(nèi)應(yīng)力對推力貢獻的標度關(guān)系及其擬合參數(shù)與雷諾數(shù)的關(guān)系Fig.13 Scaling law of the contribution of the internal stress to thrust and the relationship between fitting parameters and Reynolds number

3 結(jié)論

為了研究鲹科魚類的推進與流場結(jié)構(gòu)的關(guān)系,基于虛功率原理,將推力分解為魚體邊界變速運動的貢獻FT,a、流場中流體旋轉(zhuǎn)和應(yīng)變速率相對大小的貢獻FT,Q、壁面剪切應(yīng)力的類摩阻分量FT,v和壁面摩阻分量FT,f.除了流體旋轉(zhuǎn)和應(yīng)變速率相對大小對推力的貢獻是流場體積分外,其余3 個推力分量均為魚體表面積分,而壁面加速度對推力的貢獻是由魚體運動瞬時決定的.

通過對鲹科魚類波動的二維模型的推進分析,結(jié)果顯示推力主要來源于魚體邊界變速運動的貢獻,它是黏性流動中流體對魚體波動的瞬時反作用力,而流場中流體旋轉(zhuǎn)和應(yīng)變速率相對大小的貢獻(時均值)總是負的,它們是推力組成的主要兩個部分.從流場分區(qū)來看,推力主要來源也是魚尾的變速運動和魚尾邊界層的流動,而流場尾跡的旋渦對推力的影響較小,距離越遠影響越小.

進一步研究表明,推力分解的4 個量均與表征非定常性的Strouhal 數(shù)呈二次關(guān)系;邊界變速運動對推力的時均貢獻與雷諾數(shù)無關(guān),其余3 個量均受雷諾數(shù)影響.結(jié)合標度律分析后發(fā)現(xiàn),流場中流體的旋轉(zhuǎn)和應(yīng)變速率相對大小對推力的貢獻也有一部分是與雷諾數(shù)無關(guān),其與雷諾數(shù)有關(guān)的擺動標度部分的阻力成分相當于壁面剪切力的類摩阻分量和摩阻分量之和的一半,而類摩阻和摩阻還存在常阻力成分.

總之,基于虛功率原理的鲹科魚類波狀推進分析表明,魚體的推進力主要依賴尾部的變速運動和尾部兩側(cè)的邊界層控制,且提高波速和增強擺動速率都可以獲得高的推進力.

附錄A 數(shù)值結(jié)果驗證

本文基于OpenFOAM 通過有限體積法求解二維非定常不可壓縮Navier-Stokes 方程組(方程(3)),采用PIMPLE 算法求解控制方程.PIMPLE 算法是基于SIMPLE 算法和PISO 算法的一種改進算法,在時間步內(nèi),通過低松弛穩(wěn)定計算獲得較大的時間步長,在每個時間步內(nèi)使用PISO 算法對壓力進行多次修正使之滿足壓力速度耦合方程[B1].方程中的時間導數(shù)項采用一階隱式Euler 格式離散,動量方程中的對流項采用Gauss linear upwind 格式進行離散,擴散項采用Gauss linear 格式進行離散.

本文設(shè)定的計算域為 16L×8L,網(wǎng)格是使用OpenFOAM自帶的網(wǎng)格自動生成工具snappyHex Mesh 生成,魚體附近的動網(wǎng)格如圖A1 所示.本文計算了在 λ=1.0,c=2.0,Re=5000時,推力系數(shù)和側(cè)向力系數(shù)隨時間的變化,如圖A2 所示.網(wǎng)格選取了 160×80,480×240 以及 800×400 這3 種情況,時間步長選取了 ΔT=2.5×10-4,ΔT=1.0×10-4以及 ΔT=5.0×10-5.從圖A2 中可以清楚地看出,160×80 和 480×240 網(wǎng)格大小的計算結(jié)果存在合理的差異,而 480×240 和 800×400 網(wǎng)格大小的差異幾乎可以忽略不計.此外,計算結(jié)果與Dong 等[12]公布的結(jié)果吻合較好.因此,本文中的模擬都使用了 480×240 的中等網(wǎng)格尺寸和時間步長 ΔT=1.0× 10-4,魚體附近的最小網(wǎng)格高度約為 0.0001.

附圖 A1 魚體附近動網(wǎng)格示意圖Fig.A1 Schematic of dynamic meshes near fish body

附圖 A2 不同網(wǎng)格和時間步長的數(shù)值結(jié)果比較Fig.A2 Comparison of numerical results with different meshes and time steps

附錄B 參考文獻

B1 徐濤濤,胡文蓉.鱗鲀倒游機動的機理研究.水動力學研究與進展:A 輯,2022,37(3): 277-283 (Xu Taotao,Hu Wenrong.Study on mechanism of balistid’s backward maneuvering.Chinese Journal of Hydrodynamics,2022,37(3): 277-283 (in Chinese))