楊國黨, 胡 曉*， 石小濤， 余英俊, 張 奔, 龍澤宇

(1.三峽大學 水利與環(huán)境學院，宜昌 443000； 2.武漢大學 水利水電學院，武漢 430072)

1 引 言

自然界中魚類具有較高的機動性，與人造機器魚相比，其具有較高的推進效率。魚可以在三維空間中自由游泳，靈活躲避捕食者或者搜尋獵物。對魚類運動學的量化是分析魚類運動力學的基礎。

魚在水中游泳復雜多樣，不同水環(huán)境中呈現(xiàn)不同游泳姿態(tài)。為捕捉魚的運動，通常以高速錄像的形式獲取魚游泳行為，從中提取其運動學變量[1]。為了量化魚在生育發(fā)展各個階段的變化，F(xiàn)ontaine等[2]開發(fā)了一個視覺跟蹤系統(tǒng)來估計魚的姿勢?？律钡萚3]基于Matlab軟件開發(fā)了一種跟蹤魚類中線的程序來獲取鰱魚的擺尾頻率、加速度和速度等運動指標。為了進一步描述普遍魚類運動學規(guī)律，Liu等[4]將魚體的波動用波狀函數(shù)描述，后續(xù)有很多學者對魚波動方程進行了改進優(yōu)化。

為了實現(xiàn)有效的運動控制，魚類能夠實現(xiàn)身體變形，且由此產(chǎn)生的速度和生物力學相聯(lián)系，從而實現(xiàn)其高效游泳。魚在結構復雜的棲息地游泳需要較高的機動性，包括其方向的改變程度和完成動作所需的空間[5]。有研究發(fā)現(xiàn)，較小的魚具有較高的機動性，這主要得益于較小的魚具有較高的彎曲程度[6]。沈昊嫣等[7]探究了魚體曲率波的傳播速度，發(fā)現(xiàn)主動彎矩的驅動頻率高于魚體結構基頻時，曲率波的波速小于驅動波的波速。有研究提到，力學可以深化對魚體形態(tài)變化和運動方式之間的理解[8]。王仲威等[9]基于曲率模型發(fā)現(xiàn)斑馬魚在彎曲階段存在較大的正向力矩，而回擺階段存在比較小的反向力矩。Thandiackal等[10]研究發(fā)現(xiàn)斑馬魚在轉彎的早期階段，彎曲接近最大曲率時，魚體尾部區(qū)域壓力作負機械功。由此可見，魚類游泳姿態(tài)和力學機制對于解析魚類行為十分重要。

綜上所述，國內(nèi)外很多學者對魚類的行為學及動力學有所研究，但很少對魚自由游泳時曲率和周身壓力進行量化研究。因此本文以草魚幼魚為研究對象，基于實驗擬合振幅曲線并驗證擬合結果，量化了魚體中線曲率和輪廓曲率的變化，分析了周期內(nèi)曲率和壓力的變化趨勢和相關關系。本文可增加對鲹科魚類自游泳彎曲狀態(tài)和力學機制的認識并提供參考資料。

2 實驗設計與方法

2.1 實驗對象選取

為了獲取真實魚的運動學參數(shù)，以草魚幼魚為研究對象。實驗對象購置于湖北省宜昌市某水產(chǎn)市場，選取體質活力較好的20尾進行實驗研究。體長平均長5.2 cm±0.2 cm。實驗前將實驗對象放置于圓形玻璃水缸中暫養(yǎng)，水溫保持在18 ℃±2 ℃，實驗前24 h禁止喂食，全天不間斷向魚缸中充氧以保證溶解氧不小于6 mg/L。

2.2 實驗裝置與方法

實驗裝置如圖1所示，在透明玻璃水槽中進行。該玻璃水槽的實際尺寸為長×寬×高=34.5 cm×20 cm×22 cm。實驗前在玻璃水槽上方平行架設高速攝像機。水槽右側架設激光拍攝系統(tǒng)，激光筆采用諾青激光筆(型號為NQ-506，波長為532 nm，功率為100 mW)，棱鏡為D形柱。圖像處理系統(tǒng)選用筆記本電腦。示蹤粒子密度為1.02 g/cm3，直徑為20 μm～30 μm[11]。

圖1 實驗裝置Fig.1 Experimental device

2.3 視場獲取方法

實驗在無外界干擾的環(huán)境中進行。實驗前在水槽中均勻撒入示蹤粒子。實驗時激光束通過D形棱鏡形成約為2 mm厚的片光源，通過調(diào)節(jié)支架高度使片光源能夠照射在魚體側線上。實驗前嚴禁打開激光筆，以免對魚造成驚嚇。另外為了能夠較好地看到魚的側線變化情況，在玻璃水箱一定高度架設隔板，防止草魚幼魚在空間上做俯仰運動，另外為防止激光反射形成亮斑，在水槽側壁粘貼不透明膠帶，減少干擾。高速攝像機空間分辨率為640×480(pixels2)，每個視場拍攝5 min。

2.4 基于速度場的壓力場重構基本方法原理

采用直接積分進行壓力重構，基于二維不可壓縮Navier-Stokes方程對壓力梯度項進行積分。

不可壓縮N-S方程為

(1)

兩點之間壓力差按式(1)確定。如x1和x2兩點之間壓力差為

(2)

又因為流體連續(xù)性方程為

(3)

式中xi為坐標，ui為xi方向上的速度(i=1，2，3，分別為x，y和z方向)。

由式(3)可知，壓強梯度為粘性項和物質加速度之和，其中物質加速度歐拉表達式為

(4)

粘性項的速度場節(jié)點采用中心差分格式，表達式為

(5)

(6)

可以通過實測得到上述方程右側節(jié)點處的速度，并結合式(3)通過兩點之間的積分算出壓力值。在積分時，測量誤差是x1～x2路徑累積，當前通過沿8條路徑積分取其平均值，可大幅減小積分誤差，以實現(xiàn)準確的對壓力估算[12]。8條積分路徑如圖2所示，每條起源于域邊界，分別從左上角(LU)、上(U)、右上角(RU)、左(L)、右(R)、左下角(LD)、下(D)和右下角(RD)向每個方向傳播，平均各路徑后則為點z的壓強值。

(pL+ΔpL)+(pLU+ΔpLU)+(pLD+ΔpLD)+

(pRU+ΔpRU)+(pRD+ΔpRD)]

(7)

圖2 積分路徑Fig.2 Integral path

2.5 試驗圖像處理方法

采用互相關軟件PIVlab 2.02處理速度場，利用標準FFT相關算法，查詢區(qū)窗口大小為64×32(pixels2)。在2D視場下速度場分析之前對魚可視輪廓遮罩，后續(xù)運用基于Matlab的運行程序QUEEN2計算壓力場。為了便于后面模擬仿真，忽略胸鰭和簡化了尾鰭，對魚體主要輪廓進行分析。為了得到精細的魚體周身流場信息，參考文獻[13]研究理論，將幼魚魚體劃分為頭部、中部和尾部，具體定義為魚吻到胸鰭之間的區(qū)域為頭部(0～1/5BL)；胸鰭與腹鰭之間區(qū)域為中部(1/5～3/5BL)；腹鰭到尾柄區(qū)域為尾部(3/5～1BL)，如圖3所示。

圖3 魚體劃分Fig.3 Fish division

2.6 運動學分析

2.6.1 魚體波動振幅

魚類游泳具有較高的機動性和靈活性，部分學者將魚類的游動分為月牙尾推進、鰻鱺模式和鲹科模式[14]。草魚幼魚波動主要集中在后半部分，為鲹科模式。為了便于分析運動，此處定義草魚幼魚擺尾向前波動的過程中，當魚尾向右擺動最大幅度時再擺動一個來回，這個過程歷時為一個周期T。最終選擇魚體長BL=0.052 m的草魚幼魚近似直線游泳的視頻，將其逐幀分解并參考文獻[7,15]的理論方法提取魚體中線，并建立坐標系，如圖4(a)所示。

圖4 相對體長下魚體波動振幅曲線擬合Fig.4 Fitting of fish wave amplitude curve under relative body length

為便于分析草魚幼魚行為特征，對魚的體長采用無因次化處理。利用式a(x/BL)=a0+a1(x/BL)+a2(x/BL)2+a3(x/BL)3進行二次和三次擬合振幅曲線，如圖4(b)所示，擬合振幅曲線系數(shù)列入表1。

表1 振幅曲線擬合系數(shù)Tab.1 Amplitude curve fitting coefficients

2.6.2 魚的波動行進方程

有關研究表明，通過對魚類的二維研究，即使體型不同的魚類，如鰻魚、鱒魚、鯖魚和金槍魚等，也顯示出極其相似穩(wěn)定的波動運動[16]。鲹科模式魚類具有統(tǒng)一的運動學方程，為了形象地表示魚的運動行為，將振幅a(x/BL)、波長λ、頻率f和橫向位移y(t,x/BL)組合成參數(shù)化的數(shù)學方程。魚的運動方程可以描述為[17]

y=y0+a(x/BL)sin [2π(x/λ)-ft]

(8)

為了說明后續(xù)模擬的可操作性，利用Matlab軟件繪制相對體長下魚的包絡線,如圖5所示。圖5 的包絡線與圖4提取結果具有可比性，魚體頭部橫向位移幅度相對于尾部較小，魚的尾部具有較高的靈活性。從運動學的角度來看本文研究結果實現(xiàn)了將運動學以參數(shù)化方程來描述，這為后面的模擬仿真提供了可行性基礎。

圖5 相對體長下魚體包絡線Fig.5 Envelope diagram of fish relative to body length

2.7 計算方法

2.7.1 數(shù)值計算

使用帶有可控運動程序的水翼來復制柔性鰭或魚體是研究魚類運動學和推進力的有力工具[18,19]。為了便于仿真，本文以NACA0012作為魚的二維簡化模型，按實際魚等比例縮放建立物理模型。對整個計算流體域進行網(wǎng)格離散劃分，對魚體游泳區(qū)域進行網(wǎng)格加密，單元總數(shù)為128684，網(wǎng)格模型如圖6所示。進口采用速度進口(ν=0)，出口設置為壓力出口。速度-壓力耦合采用SIMPLE算法進行修正。為了達到較高的收斂性，時間步長取0.0005 s。鑒于數(shù)值方法在文獻[20]已有介紹，此處不再詳細闡述。

圖6 網(wǎng)格Fig.6 Grid

2.7.2 波動曲率

曲率是一種魚體內(nèi)部一致的測量，不需要完整的身體幾何標志，表現(xiàn)了魚在游泳過程中整個身體彎曲程度的度量[21]。曲率可表示為

(9)

2.7.3 魚周身壓力計算

魚擺尾游動時與流體相互作用產(chǎn)生壓力場。本文聚焦分析一個擺尾周期中線、體表曲率和魚周身壓力的相互關系。流體作用在魚體上的壓力可以表示為[22]

(10)

式中n為法向量，垂直體表指向外方向，p為壓強值，A為輪廓長度上單位厚度的面積。

3 結果與分析

3.1 模型驗證

為了進一步說明模擬結果的準確性和分析的可行性，選取t=0.65 s以及t=1.3 s時的游泳姿態(tài)，利用Tecplot與Origin對魚體周身的壓力分布進行分析，如圖7所示。結果顯示三次擬合振幅曲線比二次擬合振幅曲線能夠更好地控制幼魚擺尾行為，同時魚體的周身流體壓力與實測值的吻合度更高，因此，本文采用三次擬合的振幅曲線來模擬幼魚魚體的包絡線運動規(guī)律，在此基礎上對模擬結果進行分析。

圖7 魚周身實驗值與仿真模擬壓力值對比Fig.7 Comparison of experimental value and simulated pressure value of fish body

3.2 魚體波動曲率變化規(guī)律分析

3.2.1 不同游泳姿態(tài)下魚體中線曲率分析

圖8為魚游泳行進中線曲率變化分布。初始時刻，曲率最大處約在魚體尾部前端，對應時刻曲率最大值為kmax=29 m-1，如圖8(a)所示。從圖8(b～d)可以看出，對應時刻曲率的最大值由尾部前端向后端發(fā)展。當t/T=1/2時，魚的曲率最大處由尾部后端又逐漸發(fā)展到尾部前端。t/T=5/8～1是魚回擺階段，對應時刻kmax變化與之前相似，如圖8(f～h)所示。這表明在一個周期內(nèi)，草魚幼魚中線的kmax主要限于尾部前端和尾部后端之間往復變化。曲率最小值無明顯變化規(guī)律。

總體上，在t/T=1/2之前，魚自由前進過程中對應時刻最大曲率kmax呈先增大后減小變化趨勢；在t/T=1/2之后，對應時刻曲率最大值與之前變化相同，方向相反。這在一定程度上說明魚尾在擺動過程中草魚幼魚后半部分并非是向后伸展變化，而是向前有一定的彎曲?？傮w而言，對應時刻kmax分別在t/T=0，t/T=1/2和t/T=1時相對達到最小，而在t/T=1/4和t/T=3/4時達到相對總體最大值。當t/T=1/4～3/4，對應時刻kmax呈現(xiàn)先減小后增大的原因是魚正處于S型擺尾階段[23]，中線彎曲最大處主要集中在魚體后部。整個周期中對應時刻曲率最大值相對總體最大值和最小值分別為102 m-1和29 m-1。

有研究提到魚體曲率與角度變化有關，但是很少提到魚體變形與力學關聯(lián)機制[24]。為了說明魚體變形影響機制，本文將魚體按體長劃分為 10段。圖9為每段中間點最大曲率在周期內(nèi)的變化趨勢，魚體彎曲最大處主要集中在魚體尾部，這為3.3節(jié)壓力分析奠定基礎。

圖8 不同時刻魚體中線曲率變化Fig.8 Curvature changes of fish body at different times

圖9 周期內(nèi)魚體每段中線曲率峰值Fig.9 Peak curvature of each central line of the fish body in a period

3.2.2 中線曲率與周身曲率耦合

魚頭部分較硬，可視作剛體，魚身部分富有彈性，可視作彈性體[25]。因為魚體頭部中線曲率和對應體表曲率的不一致性，本文只分析如圖10所示魚體中部和尾部對應中線曲率與周身曲率的關系。在t/T=1/8～1/2時，分析發(fā)現(xiàn)除魚頭部分外，不同游泳姿態(tài)下魚體中線的曲率和幼魚周身曲率均符合K=ak+b線性函數(shù)關系(k為中線曲率，K為周身輪廓曲率，P<0.05)。

圖10 幼魚中部和尾部的中線曲率與體表輪廓曲率關系Fig.10 Young fish line at the end of the central,curvature relationship with surface contour curvature

3.3 幼魚周身流體壓力

3.3.1 魚周身壓力分布

圖11為幼魚周身的流體壓力云圖。如圖11(a)所示，初始時刻魚頭周圍流體表現(xiàn)為負壓，魚尾下側表現(xiàn)正壓，上側表現(xiàn)負壓，此時合壓力變化范圍為-9.6 Pa/BL～0 Pa/BL。之后魚頭向下側輕微擺動，表現(xiàn)為下側流場為正壓，上側表現(xiàn)負壓且有增大趨勢，如圖11(b)所示。t/T=1/4時，如 圖11(c)所示，魚體頭部有大面積高壓力區(qū)出現(xiàn)，這是由魚頭向下擺動速度誘導的，此刻魚尾壓差變化區(qū)間為-1.9 Pa/BL～0 Pa/BL。之后在慣性力作用下魚尾繼續(xù)向下擺動，顯而易見在t/T=3/8時魚尾表面流體壓力發(fā)生變化，魚尾下側的流體壓力高于右側。當t/T=1/2時，在幼魚自身行為主導驅動下，魚體下側的擺尾幅度達到極限，如圖11(e)所示。t/T=5/8～1時，魚尾開始回擺，整個壓力場變化與前半個周期相反，如圖11(f～i)所示。

圖11 魚周身壓力變化云圖Fig.11 Cloud chart of pressure changes around fish

值得注意的是，在t/T=1/4和t/T=3/4時刻，約在魚體尾部流體表現(xiàn)為正負壓交替現(xiàn)象，主要原因是各部位的擺動速度和彎曲程度不同所致。

3.3.2 魚體分段曲率和壓力變化分析

為了探尋魚體每段中線曲率與對應魚體表面的流體壓力關系，擬合每一段中線曲率和魚周身流體壓力，如圖12所示。從圖12(a～b)可以看出，魚體頭部每段中線曲率與下側表面流體壓力變化趨勢基本一致；相應地在魚體尾部，每段中線曲率主要與魚體上側表面流體壓力隨時間變化趨勢一致(圖12(g～i))；然而在魚體中部約(1/2～3/5)BL區(qū)間段，中線曲率和魚體上下表面流體壓力均未隨時間呈現(xiàn)一致的變化趨勢(圖12(f))。此外，魚體表面流體壓力與曲率在周期內(nèi)均出現(xiàn)交叉變化，曲率波峰最大滯后壓力波峰約為3/8T。結合圖9可知，每段中線曲率最大值變幅越大，周期內(nèi)相位滯后隨魚體長度的變化就越明顯，這與曲率波滯后機電信號類似[26]。

圖12 周期內(nèi)分段曲率和壓力隨時間變化曲線Fig.12 Curve of piecewise curvature and pressure over time in a period

3.3.3 魚周身壓力-曲率耦合

圖13顯示了魚體中部和尾部的表面流體壓力與曲率的關系。由于t/T=5/8～1是魚回擺過程，因此只選擇t/T=1/8～1/2區(qū)間作為展示。從圖13(a)可以看出，魚體中部流體壓力主要隨曲率增大而增大(P<0.05)。隨著魚體向前游動，魚體中部和尾部的上下側流體壓力均隨曲率的增大而增大(P<0.05)，如圖13(b)所示。當t/T=3/8時，魚體中部下側流體壓力和曲率呈正相關關系，而上側流體壓力與曲率呈負相關變化關系，如 圖13(c)所示；魚體尾部上下側流體壓力和曲率均呈正相關關系(P<0.05)。當魚體向下擺尾最大限度時，魚體中部兩側表面壓力和曲率呈負相關關系，魚體尾部表面流體壓力和曲率表現(xiàn)為負相關關系(P<0.05)，如圖13(d)所示。

從圖13可以看出，魚體中部和尾部兩側表面流體壓力和曲率大多表現(xiàn)出明顯的線性關系(K=aF+b)。結合圖11可知，流體與魚體相互作用時，魚對流體的壓力像波一樣向四周擴散，并逐漸減弱，所以魚對流體壓力并非持續(xù)性增大，又因為曲率最大處并非是壓力最大處，所以會出現(xiàn)曲率和壓力呈負相關變化。另外由于魚逐漸擺尾向下游動，魚體中部和尾部之間出現(xiàn)負壓和正壓逐漸過渡，這也是曲率較大處呈現(xiàn)壓力先減小后增大(或先增大后減小)的原因。此外，由于幼魚中間區(qū)域和尾部中線曲率與表面流體壓力有較明顯的線性關系(P<0.05)，也說明魚體中部和尾部整體表面流體壓力與中線曲率也呈線性關系。

圖13 魚體中部和尾部周身流體壓力與魚體曲率的關系Fig.13 Correlation between pressure and curvature in the middle and tail of fish

4 結 論

本文利用仿真模擬的方法分析了草魚幼魚的自由擺尾行為，量化了魚體中線曲率在一個擺尾周期內(nèi)的變化，同時對幼魚周身流體壓力與魚體表面和中線曲率的相關性進行了分析，探究了幼魚中部和尾部的周身流體壓力和曲率的變化規(guī)律，主要結論如下。

(1) 擬合三次振幅曲線方程比傳統(tǒng)二次振幅曲線方程能夠更好控制幼魚的自由擺動，且得到的魚體周身流體壓力與PIV實測結果吻合度更高，相比常用的二次模擬結果要更精確。

(2) 幼魚中部和尾部的中線曲率與體表曲率呈較為明顯的線性關系(P<0.05)，且魚體中線最大曲率的變化主要局限于幼魚尾部。

(3) 幼魚尾部和頭部的中線曲率分別與魚體上下側的流體壓力變化趨勢一致；而在幼魚的中間部位，魚體的中線曲率與流體壓力未有一致走勢，曲率波峰滯后壓力波峰最大可達3/8T。

(4) 幼魚中部和尾部的周身流體壓力在前半個運動周期內(nèi)和魚體體表、中線的曲率均呈線性相關。

本文限制了幼魚垂向運動，通過實驗獲取幼魚運動學參數(shù)和周身截面流場。利用實驗數(shù)據(jù)擬合分析及標定數(shù)值模型，實驗與模擬結果的比對說明本方法可較好地模擬魚體曲率變化和周身壓力的關系，對理解幼魚自推進變形和受力機制具有一定的意義。