尾鰭擺動推進水動力與流動渦結(jié)構(gòu)特征分析

蘇博越，翟樹成，白亞強，張 軍

（中國船舶科學研究中心，江蘇 無錫 214082）

0 引 言

大多數(shù)魚類通過身體和（或）尾鰭的擺動，推動水流，使魚體向前運動，人們將這種運動仿生定義為BCF 模式（body-caudal fin）。這種模式通常具有游動速度快、推進效率高的特點，從而成為仿生推進研究的熱點。1994年，MIT研究制造出世界上第一條機器魚，仿金槍魚機器魚Robotuna[1]，此后各國高校與科研院所紛紛開展仿生機器魚研制和仿生推進技術(shù)研究。英國ESSEX 大學的胡豁生團隊[2-6]研制出了G 系列機器魚，能以90°/0.2 s 的最大速度實現(xiàn)大角度、高頻率的運動。2009 年，哈爾濱工程大學[7]以藍鰭金槍魚為藍本研制出仿生水下機器魚“仿生-Ⅰ號”，依靠月牙形尾鰭擺動獲取推力，尾鰭擺動頻率為1.33 Hz，游動速度最高可達1.2 m/s。隨著計算機科學、計算方法和流體仿真軟件的發(fā)展，數(shù)值模擬和仿真方法已成為研究魚類復雜游動現(xiàn)象的一個有效手段。日本Chiba大學的劉浩[8-10]在20世紀90 年代中后期引入人工偽壓縮因子，采用有限體積法對Navier-Stokes 方程進行離散化，實現(xiàn)計算流體動力學（CFD）的解算、流場顯示及動力分析等后處理操作，并將CFD 方法成功應用于幼魚游動、飛蛾懸停等生物游動、飛行的復雜流動機理研究。Wolfgang[11]等結(jié)合活魚的試驗數(shù)據(jù)，使用三維面元法計算了金槍魚的游動，并分析了尾流中渦旋之間的相互作用。夏全新（2005）[12]采用有限體積法和動網(wǎng)格技術(shù)，對魚類的波動游動方式進行數(shù)值模擬。通過對二維零厚度板的波狀擺動推進進行計算，研究了在高雷諾數(shù)魚類勻速運動狀態(tài)下的流場分布，并在尾跡中發(fā)現(xiàn)了“之”字形排列的渦環(huán)結(jié)構(gòu)。2012 年張曦[13]對藍鰭金槍魚尾鰭擺動推進水動力性能進行了多角度的分析，包括在振動半圓柱擾流場中二維振蕩水翼的推進性能、尾鰭形狀和尾鰭柔性以及非對稱運動對擺尾推進性能的影響等。2015年，白亞強等[14]基于彈性光順模型和局部網(wǎng)格重構(gòu)模型動網(wǎng)格技術(shù)，建立了在大擺幅運動下柔性長鰭扭波推進水動力及流場渦結(jié)構(gòu)的數(shù)值計算方法，并通過水動力試驗從流場和水動力方面驗證了數(shù)值計算方法的正確性。2017 年李寧宇等[15]采用基于RANS 方程的數(shù)值模擬方法研究了金槍魚自主游動的水動力特性，并通過模擬金槍魚運動姿態(tài)以及巡游狀態(tài)受到的阻力與游動速度的關(guān)系，求解了金槍魚從靜止到巡游狀態(tài)下的自主游動過程。以上研究大多著眼于宏觀水動力結(jié)果，對魚游流場的精細化空間渦結(jié)構(gòu)的計算和渦結(jié)構(gòu)與水動力相互關(guān)系的研究尚不夠深入。2021年張棟[16]對近年生物游動和飛行等過程中的渦結(jié)構(gòu)模擬及力學特征規(guī)律進行了綜述，并提出了開展真實雷諾數(shù)下生物游動飛行精細數(shù)值模擬和渦結(jié)構(gòu)對生物體游動與飛行過程中受力貢獻方面研究的展望。

本文主要從BCF 推進運動模式出發(fā)，研究建立仿金槍魚尾鰭擺動推進粘性流場精細化數(shù)值計算方法以及水動力、渦結(jié)構(gòu)提取與分析方法，并通過網(wǎng)格無關(guān)性驗證及與相關(guān)文獻結(jié)果的比較等驗證本文方法的有效性。采用本方法數(shù)值研究不同運動參數(shù)組合下尾鰭擺動推進水動力隨St數(shù)的變化規(guī)律，結(jié)合尾鰭擺動運動、推力、表面壓力等，分析尾鰭擺動流動渦結(jié)構(gòu)特征及其演化過程。

1 計算方法介紹

1.1 BCF運動模式

在BCF 推進模式中，尾鰭運動通常由圍繞尾柄末端的搖擺和由尾柄驅(qū)動的橫向振蕩運動組成。計算模型的坐標系如圖1所示。來流速度的方向與x軸的正方向相同。尾鰭端部沿z方向橫向往復移動，同時圍繞通過尾鰭端部平行于y軸的旋轉(zhuǎn)軸線搖擺[17]。

圖1 尾鰭運動坐標系Fig.1 Coordinate system of caudal fin motions

尾鰭橫向運動的位移z(t)和尾鰭運動坐標系中尾鰭軸線的搖擺角度θ(t)滿足正弦變化規(guī)律，即

式中：A為橫移運動的幅值；θ0為搖擺運動的角幅度；f為尾鰭擺動頻率，在尾鰭擺動過程中，橫移運動和搖擺運動的頻率相等；?0為橫移運動和搖擺運動的相位差。

在尾鰭運動過程中，影響擺動尾鰭水動力性能的無量綱參數(shù)——斯特羅哈爾數(shù)（St）表示為

式中：A0表示尾鰭運動過程中脫落尾渦的寬度，近似表達為橫移運動幅值的兩倍，即A0=2A；V0為來流速度。

通過數(shù)值計算可以得到尾鰭在不同運動狀態(tài)下的水動力系數(shù)。水動力系數(shù)表示為

式中：Fx、Fz分別為尾鰭所受到的x方向、z方向的分力，F(xiàn)x方向與推力Ft方向相反，將推力系數(shù)記為Kt，則Kt=-Cx；S為尾鰭弦長與展長的乘積；C0為尾鰭弦長。

1.2 數(shù)值計算

本文采用商用軟件Fluent V17.0進行數(shù)值模擬，應用動網(wǎng)格技術(shù)來處理尾鰭運動邊界隨時間運動變形的問題。本文尾鰭邊界的運動由UDF程序控制。數(shù)值計算中網(wǎng)格動態(tài)變化過程采用彈簧光順模型（spring-based smoothing）和局部重構(gòu)模型（local remeshing）。

1.2.1 控制方程

本文采用不可壓縮粘性流體的雷諾平均方程（RANS）作為流體控制方程，方程如下：連續(xù)性方程為

動量方程為

式中，ui為平均速度分量（u1=u,u2=v,u3=w），ρ為流體質(zhì)量密度，p為平均壓力，μ為流體分子粘性系數(shù)為雷諾應力項。

湍流模型采用SSTk-ω兩方程模型[18]。

本文對流項和空間項均采用二階精度格式進行離散。

1.2.2 彈簧光順模型

彈簧光順模型基于彈簧的平滑，將網(wǎng)格中任意兩個連接的網(wǎng)格節(jié)點之間的邊看作是連接兩點的彈簧，所有網(wǎng)格節(jié)點的邊形成彈簧網(wǎng)絡(luò)。邊界運動前，各個節(jié)點的初始間距為網(wǎng)格的平衡狀態(tài)，當某一邊界節(jié)點發(fā)生位移時，系統(tǒng)將產(chǎn)生一個與該節(jié)點連接的所有彈簧位移成正比的力，會引起整個彈簧網(wǎng)絡(luò)的變化，通過求解彈簧網(wǎng)絡(luò)重新達到平衡時節(jié)點的新位置即可相應地得到新的網(wǎng)格節(jié)點位置。

在尾鰭運動過程中，由于運動范圍較大，運動速度快，尾鰭表面的網(wǎng)格很容易產(chǎn)生嚴重的變形。在數(shù)值模擬過程中，通過對彈簧剛度系數(shù)的不斷調(diào)整，可以把尾鰭周圍的網(wǎng)格變形傳遞到離尾鰭較遠的網(wǎng)格，從而避免尾鰭周圍網(wǎng)格的較大變形，保證周圍網(wǎng)格的質(zhì)量。

1.2.3 局部網(wǎng)格重構(gòu)模型

彈簧光順模型通常用于網(wǎng)格區(qū)域中的動態(tài)網(wǎng)格生成。但如果運動邊界的位移遠大于網(wǎng)格尺寸，或者邊界運動方向造成網(wǎng)格的扭曲，即使采用彈簧光順模型仍可能導致網(wǎng)格質(zhì)量下降，甚至出現(xiàn)負體積網(wǎng)格，導致網(wǎng)格失真過大引起計算不收斂。采用局部網(wǎng)格重構(gòu)技術(shù)可以很好地控制網(wǎng)格質(zhì)量，防止因位移變化過大而產(chǎn)生負體積網(wǎng)格。

1.3 計算參數(shù)選取與邊界條件設(shè)置

本文計算網(wǎng)格如圖2 所示，在尾鰭附近建立加密區(qū)域（黃色區(qū)域），其尺寸為5C0×5C0×7.5C0，外圍流場尺寸為25C0×12.5C0×12.5C0。圖中藍色表面為速度入口，紅色表面為壓力出口，加密區(qū)域內(nèi)尾鰭壁面設(shè)置為無滑移邊界。

圖2 網(wǎng)格劃分Fig.2 Mesh distributions

尾鰭擺動推進流場為非定常流動，選取時間步長為0.005 s，每個時間步最大迭代次數(shù)為50，總時間步數(shù)為4000。

1.4 網(wǎng)格無關(guān)性驗證

在網(wǎng)格無關(guān)性驗證的數(shù)值計算中，尾鰭型值選自文獻[13]，其中尾鰭弦長C0=129 mm，展長D=241 mm，尾鰭橫移幅值A(chǔ)=100 mm，搖擺角幅度θ0=25°，橫移運動和搖擺運動的相位差?0=-π/2，運動頻率f=0.5 Hz，來流速度V0=0.31 m/s。

本文首先設(shè)計了網(wǎng)格稀疏程度不同的三套網(wǎng)格G1、G2、G3，其網(wǎng)格數(shù)分別為120 萬、286 萬和541萬，其加密區(qū)域網(wǎng)格數(shù)分別為90萬、137萬和185萬。通過對比不同網(wǎng)格數(shù)下的計算結(jié)果（見表1、圖3）可知，采用網(wǎng)格G1的計算結(jié)果與其它兩個結(jié)果相差較大，而采用網(wǎng)格G3的水動力系數(shù)均值計算結(jié)果與網(wǎng)格數(shù)量相關(guān)性較小。其中網(wǎng)格G3 的Cxm計算結(jié)果與網(wǎng)格G2 的結(jié)果偏差小于4%，而對于側(cè)向力系數(shù)均值Czm，由于尾鰭運動的對稱性，在數(shù)值計算達到穩(wěn)定后，一個周期內(nèi)的Czm理論值為0，從計算結(jié)果來看，可以認為其已達到無關(guān)性要求。圖3 為三套網(wǎng)格的水動力系數(shù)達到穩(wěn)定后一個周期內(nèi)水動力系數(shù)隨時間的變化曲線。因此，后續(xù)將選取網(wǎng)格G3的水動力計算結(jié)果進行比較分析。

表1 不同數(shù)量網(wǎng)格計算結(jié)果Tab.1 Results of different meshes

圖3 不同網(wǎng)格數(shù)下水動力系數(shù)隨時間變化Fig.3 Hydrodynamic coefficients varying with time for different meshes

1.5 渦結(jié)構(gòu)提取分析方法

本文采用Q準則判據(jù)提取尾鰭擺動推進過程中流場的渦特征結(jié)構(gòu)[19]。圖4為采用網(wǎng)格G3計算的流場渦結(jié)構(gòu)Q=0.01 s-1結(jié)果?？梢钥吹?，尾流中渦環(huán)的輪廓清晰光順，較小尺度渦的特征也有保留。雖然繼續(xù)增加網(wǎng)格會使渦結(jié)構(gòu)更加精細，但同時對計算資源的要求更大，因此，本文選用G3網(wǎng)格進行尾鰭擺動推進流場數(shù)值計算，并對渦結(jié)構(gòu)進行提取分析。

圖4 尾鰭擺動推進渦結(jié)構(gòu)（網(wǎng)格G3）Fig.4 Vortex structure of caudal fin(Mesh G3)

1.6 水動力性能計算結(jié)果驗證

2012 年張曦[13]采用RANS 方法計算尾鰭擺動推進水動力，同時開展了尾鰭擺動推進水動力試驗，水動力數(shù)值計算結(jié)果得到了試驗結(jié)果的驗證。從圖5中可以看到，本文Cx數(shù)值計算結(jié)果與文獻[13]計算結(jié)果趨勢一致，量值基本吻合，峰值位置稍有差異，可見本文方法水動力的數(shù)值計算結(jié)果是合理有效的。

圖5 計算結(jié)果與文獻結(jié)果比較Fig.5 Comparison of calculation results with the results in literature

2 水動力隨St數(shù)的變化

采用前文的方法對黃鰭金槍魚尾鰭擺動推進水動力與流場進行了計算分析，黃鰭金槍魚尾鰭弦長C0為0.4 m，展長D為0.66 m（如圖1）。

為了研究St數(shù)對水動力系數(shù)的影響，本文通過系列變化來流速度、擺動頻率、橫移幅度，來實現(xiàn)St數(shù)的改變，具體計算工況參數(shù)見表2。如圖6 所示，在均勻來流下單尾鰭擺動推進水動力與流動渦結(jié)構(gòu)由參數(shù)St主控，在所研究的參數(shù)范圍內(nèi)，隨著St增大，水動力系數(shù)均值Cxm逐漸減小，也即平均推力系數(shù)Kt逐漸增大。

表2 計算工況參數(shù)Tab.2 Parameters for calculation

圖6 Cxm隨St變化曲線Fig.6 Calculated results of Cxm with different St

本文將尾鰭尾柄位置橫移為0且下一時刻向z軸正向運動的時刻作為一個周期起始，記為t=0。

圖7為計算工況3下（St=0.2）一個周期內(nèi)尾鰭水動力系數(shù)Cx、Cz、Cm瞬時變化曲線?？梢园l(fā)現(xiàn)，在t=(1/8)T、(5/8)T時刻附近，尾鰭的瞬時推力最大；在t=(3/8)T、(7/8)T時刻附近，尾鰭瞬時推力最小。由于尾鰭運動是關(guān)于x軸鏡像對稱，因而推力系數(shù)Cx變化周期為尾鰭運動周期的一半；而尾鰭運動關(guān)于z軸不存在對稱性，所以側(cè)向力系數(shù)Cz和力矩系數(shù)Cm變化周期與尾鰭運動周期相同。

圖7 一個周期內(nèi)尾鰭瞬時水動力變化曲線（工況3）Fig.7 Hydrodynamic coeffcients of caudal fin within one period

3 尾渦結(jié)構(gòu)特征

圖8 為工況3 中尾鰭擺動推進一個周期內(nèi)的渦結(jié)構(gòu)（Q=0.01 s-1）演化過程。在t=(0～2/8)T尾鰭時刻內(nèi)，尾柄向紙面外（z軸正向）運動，然后在2/8T時刻反向，在t=(2/8～6/8)T尾柄持續(xù)向紙面內(nèi)運動至位移最大處，在t=(6/8)T后改變方向向初始位置運動。由于尾柄平移和尾鰭轉(zhuǎn)動兩種運動間存在相位差，本次計算結(jié)果中尾渦的脫落時刻也相較于尾柄運動有些滯后。比較圖7 和圖8 可以發(fā)現(xiàn)，在t=(1/8)T、(5/8)T附近尾鰭水動力系數(shù)Cx、Cz達到極值，此刻尾鰭處于梢部脫落連續(xù)梢渦的階段，為渦環(huán)形成的中期；在t=(3/8)T、(7/8)T附近尾鰭水動力系數(shù)Cx、Cz為0，此時尾鰭表面渦完成脫落，并與梢渦組合形成渦環(huán)結(jié)構(gòu)，本文定義此渦環(huán)結(jié)構(gòu)為初始渦環(huán)。一個周期內(nèi)尾鰭往復運動產(chǎn)生了兩個完整的渦環(huán)，每個渦環(huán)的產(chǎn)生過程中Cx完成一個周期的變化，而由于前后半個周期渦環(huán)z向符號發(fā)生變化，Cz僅變化了半個周期。

圖8 一個周期內(nèi)尾鰭渦結(jié)構(gòu)演化（工況3）Fig.8 Vortex structure development of caudal fin within one period

圖9～12 為不同擺動頻率流場渦結(jié)構(gòu)特征（Q=0.01 s-1）隨St數(shù)變化結(jié)果。如圖10 所示，在St=0.2時，從尾鰭脫落的渦環(huán)在流場中交錯排列，靠近尾鰭的渦環(huán)尺寸較小，且相鄰渦環(huán)相互連接。隨著尾流演化，相鄰渦環(huán)開始發(fā)生分離。通過比較可以發(fā)現(xiàn)在St數(shù)較小時，渦環(huán)之間的分離較慢，渦結(jié)構(gòu)呈“之”字型；在St數(shù)較大時（St≥0.3），相鄰渦環(huán)之間的分離較快，尾流渦結(jié)構(gòu)變?yōu)殡p列楔形渦；在St數(shù)達到0.5 后，由于渦環(huán)間距縮短，雙列渦的相鄰渦環(huán)也逐漸開始相互連接，且在尾鰭梢部后方出現(xiàn)上下兩條小尺度的“之”字形帶狀渦結(jié)構(gòu)，其渦強度也隨St數(shù)增加而增大。將t=0、(4/8)T時刻，上半個周期產(chǎn)生的初始渦環(huán)與x軸正向夾角定義為初始渦環(huán)角，初始渦環(huán)角隨St數(shù)變化如圖13 所示。結(jié)合圖6、圖13 可以發(fā)現(xiàn)，在參數(shù)范圍內(nèi)，隨著St數(shù)增加初始渦環(huán)角逐漸增大，尾鰭運動產(chǎn)生的推力也逐漸增大。

圖9 擺動頻率0.5 Hz渦結(jié)構(gòu)（St=0.2）Fig.9 Vortex structure result of 0.5 Hz oscillating frequency(St=0.2)

圖10 擺動頻率0.75 Hz渦結(jié)構(gòu)（St=0.3）Fig.10 Vortex structure result of 0.75 Hz oscillating frequency(St=0.3)

圖11 擺動頻率1Hz渦結(jié)構(gòu)（St=0.4）Fig.11 Vortex structure result of 1 Hz oscillating frequency(St=0.4)

圖12 擺動頻率1.25 Hz渦結(jié)構(gòu)（St=0.5）Fig.12 Vortex structure result of 1.25Hz oscillating frequency(St=0.5)

圖13 初始渦環(huán)角隨St數(shù)變化Fig.13 Initial vortex ring angle changing with St

4 結(jié) 論

本文采用動網(wǎng)格技術(shù)，建立了均勻來流下尾鰭擺動推進流場與水動力的數(shù)值計算方法，并通過網(wǎng)格無關(guān)性計算和其它文獻結(jié)果比較驗證了計算方法。采用Q準則判據(jù)提取尾鰭擺動推進流場的渦結(jié)構(gòu)，并通過比較水動力與渦結(jié)構(gòu)的周期性變化，說明了尾鰭擺動推進水動力與渦結(jié)構(gòu)的關(guān)系。同時，本文還分析了渦結(jié)構(gòu)特征與St數(shù)的關(guān)系。研究表明：

（1）本文數(shù)值計算方法可模擬尾鰭擺動推進流場，獲取尾鰭擺動推進水動力和渦結(jié)構(gòu)的演化。尾鰭擺動推進尾流精細渦結(jié)構(gòu)，其整體上呈交錯分布的環(huán)狀結(jié)構(gòu)，并隨尾流演化逐漸分離。

（2）在均勻來流下，尾鰭擺動推進水動力由參數(shù)St主控，在所研究的參數(shù)范圍內(nèi)，隨著St增大，平均推力系數(shù)Kt逐漸增加。

（3）從一個周期渦結(jié)構(gòu)的演化來看，在梢渦持續(xù)脫落時尾鰭水動力系數(shù)Cx、Cz達到極值，在尾鰭表面渦完全脫落，尾流中出現(xiàn)完整渦環(huán)結(jié)構(gòu)時，尾鰭水動力系數(shù)Cx、Cz為0。

（4）隨著St數(shù)增加，初始渦環(huán)結(jié)構(gòu)與x軸正向夾角逐漸增大，在St=0.3左右時，尾鰭尾流中渦環(huán)結(jié)構(gòu)逐漸從單列演變?yōu)殡p列。