基于虛功率原理的鲹科魚類波狀推進分析1)

黃愷俊 余永亮

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

引言

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

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

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

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

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

1 計算模型與分析方法

1.1 魚游模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 結果與討論

2.1 典型問題的推力分解

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

作為典型的案例,對比了一個擺動周期里推力系數隨時間的變化,如圖3(a)所示.結果表明基于虛功率原理的推力曲線與由式(4)直接積分的推力曲線重合,具有一致性.進一步地,圖3(b)給出了4 個分量CT,a,CT,Q,CT,v和CT,f隨著時間變化的曲線,發現邊界加速度對推力始終是一個正的貢獻,而流場中Q對推力有很大的負貢獻,兩者相位幾乎相反.邊界的類摩阻分量相對于其他分量要小得多,幾乎可以忽略不計,而摩阻分量在一個周期內始終是阻力.四者的周期均值分別為 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 推力系數在一個周期內的變化Fig.3 The time-dependent thrust coefficient in an undulation cycle

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

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

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

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

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

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

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

圖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

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

圖6 流場分區示意圖Fig.6 Schematic of flow field zoning

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

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

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

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

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

2.2 推力與St 數關系

圖9 不同雷諾數Re條件下,〈CT,a〉,〈CT,Q〉,〈CT,v〉和〈CT,f〉 隨運動參數St的變化規律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情況下的各項貢獻的堆疊柱狀圖.從圖中可以看出,無論是低雷諾數還是高雷諾數情況,正推力的貢獻都主要來自于邊界變速運動.壁面切應力的類摩阻分量無論產生正貢獻還是負貢獻,相對于其他3 項總是較小的.阻力分量主要來源于流場中Q的貢獻和壁面摩阻的貢獻.對于低雷諾數(Re=500)的情況,兩者的阻力貢獻相當,摩阻的阻力貢獻略大于流場中Q的貢獻.對于高雷諾數(Re=5000) 的情況,流體中Q分布的阻力貢獻強于摩阻的貢獻.

圖10 不同雷諾數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)給出的分區貢獻占比不會隨著運動參數S t的變化而變化.

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

2.3 波狀推進的標度律分析

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

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

圖13 內應力對推力貢獻的標度關系及其擬合參數與雷諾數的關系Fig.13 Scaling law of the contribution of the internal stress to thrust and the relationship between fitting parameters and Reynolds number

3 結論

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

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

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

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

附錄A 數值結果驗證

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

本文設定的計算域為 16L×8L,網格是使用OpenFOAM自帶的網格自動生成工具snappyHex Mesh 生成,魚體附近的動網格如圖A1 所示.本文計算了在 λ=1.0,c=2.0,Re=5000時,推力系數和側向力系數隨時間的變化,如圖A2 所示.網格選取了 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 網格大小的計算結果存在合理的差異,而 480×240 和 800×400 網格大小的差異幾乎可以忽略不計.此外,計算結果與Dong 等[12]公布的結果吻合較好.因此,本文中的模擬都使用了 480×240 的中等網格尺寸和時間步長 ΔT=1.0× 10-4,魚體附近的最小網格高度約為 0.0001.

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

附圖 A2 不同網格和時間步長的數值結果比較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))