馮億坤，蘇玉民，徐小軍，劉煥興，王兆立

（1.國防科技大學(xué)智能科學(xué)學(xué)院，長沙 410073；2.哈爾濱工程大學(xué)水下機器人技術(shù)重點實驗室，哈爾濱 150001；3.北京特種機械研究所，北京 100143；4.北京控制工程研究所，北京 100094）

0 引 言

海洋覆蓋了地球表面的71%，蘊藏著豐富的自然資源，但海洋的平均深度達(dá)到三千多米，最深處超過一萬米，這是人類所無法直接到達(dá)的。因此，各種各樣的水下運載器應(yīng)運而生，并且得到了快速的發(fā)展。槳-舵推進(jìn)系統(tǒng)作為最成熟的水下推進(jìn)與操縱裝置，被廣泛地應(yīng)用于水下運載器的研制，然而其存在體積大、噪音高、推進(jìn)效率低以及操縱性能差的缺點。受水中魚類游動的啟發(fā)，研制具備魚類優(yōu)異游動性能的仿生水下機器魚在近三十年得到了廣泛的研究，從麻省理工學(xué)院在1994年研制成功世界上第一條真正意義上的仿生機器魚Robo-Tuna[1-2]開始，國內(nèi)外科研工作者研制了各類仿生機器魚，進(jìn)行仿生機器魚研究的首要任務(wù)是對魚類游動機理的研究。由于魚類游動性能的復(fù)雜性，采用活魚進(jìn)行試驗無法精確地得到魚體的水動力性能和流場信息，因而難以對魚類的游動機理進(jìn)行詳細(xì)的參數(shù)化研究。隨著計算流體流體力學(xué)（CFD）方法的發(fā)展，研究者對魚類的游動開展了大量的數(shù)值模擬研究。Yang 等[3]對仿生金槍魚尾鰭的水動力性能進(jìn)行了數(shù)值計算，對尾鰭的剛性和柔性進(jìn)行了對比；張曦等[4]數(shù)值計算了仿生金槍魚尾鰭、仿海豚尾鰭和仿生鯨魚尾鰭的水動力性能，并進(jìn)行了對比分析；Li和Su[5]對仿生胸鰭的推進(jìn)性能進(jìn)行了數(shù)值計算，并對胸鰭運動過程中的三維渦結(jié)構(gòu)與胸鰭表面壓力分布進(jìn)行了詳細(xì)分析；張紀(jì)華等[6]采用實驗的方法，對柔性胸鰭的水動力性能進(jìn)行了水池試驗；Zhu 等[7]采用邊界元法（BEM）對仿生魚在均勻來流中的水動力性能進(jìn)行了數(shù)值計算；Yang 和Su[8]對仿生金槍魚在均勻來流中水動力性能進(jìn)行了數(shù)值計算，并分析了流場渦結(jié)構(gòu)?？紤]到自然界中的魚類是在自身推進(jìn)下的自主游動，Borazjani 和Sotiropoulos[9]對仿生魚的單自由度（1-DoF）自主游動進(jìn)行了數(shù)值計算。由于魚類在進(jìn)行周期性游動時，身體軀干及尾鰭的運動是周期對稱的，因此可以看作是尾鰭運動平面內(nèi)的自主游動；Xia等[10]、Li等[11]和Feng等[12]對仿生金槍魚在水平面內(nèi)的三自由度（3-DoF）自主游動進(jìn)行了數(shù)值計算，并且分析了仿生魚產(chǎn)生的三維尾渦結(jié)構(gòu)以及尾鰭表面的前緣渦（LEV）在推力產(chǎn)生中的作用；劉煥興等[13]對仿生機器魚的擺動-滑行游動模式進(jìn)行了數(shù)值計算，結(jié)果表明滑行比大于等于0.6時，擺動-滑行游動模式比連續(xù)擺尾游動模式消耗更少的能量。

考慮到魚類游動的復(fù)雜性，通常為多個方向的平移/旋轉(zhuǎn)運動耦合而成的剛性運動或者三維空間中的柔性運動，因而難以針對CFD 方法在求解魚類游動性能的有效性開展實驗驗證。本文首先建立適用于求解仿生水動力學(xué)問題的CFD 網(wǎng)格劃分策略，其次基于CFD 方法建立仿生魚水動力學(xué)中的關(guān)鍵運動部件和游動問題的數(shù)值計算方法，最后針對仿生尾鰭、仿生胸鰭和仿生魚在均勻來流中的水動力性能以及仿生魚的3-DoF 自主游動性能展開數(shù)值計算，并與文獻(xiàn)中的實驗與數(shù)值計算結(jié)果進(jìn)行對比，驗證計算方法的有效性。本文系統(tǒng)地給出了基于CFD 二次開發(fā)的仿生水動力學(xué)數(shù)值計算方法以及數(shù)值算例和驗證結(jié)果，對仿生水動力學(xué)在工程中的應(yīng)用具有重要的意義。

1 CFD數(shù)值計算方法

1.1 控制方程

自然界中魚類游動的介質(zhì)為不可壓縮粘性流體，其運動滿足時均連續(xù)性方程和雷諾平均納維-斯克托斯(RANS)方程：

式中，ui和uj為流體速度分量的時均值（i，j=1，2，3），p為流體壓力時均值，μ為流體動黏性系數(shù)，ui'uj'為速度脈動分量乘積的時均值，Si為廣義源項。

采用通用CFD 軟件包FLUENT 對流場進(jìn)行求解，基于有限體積法對式（1）和式（2）進(jìn)行離散，湍流模型采用SSTk-ω兩方程模型，選擇基于壓力的求解器，采用標(biāo)準(zhǔn)壁面函數(shù)法對近壁面流動進(jìn)行處理，采用SIMPLE 算法對壓力和速度進(jìn)行耦合，離散方式中的梯度項選擇Least Squares Cell Based，壓力項選擇Second Order，動量項選擇Second Order Upwind，時間項采用1st-Order Implict，計算物體表面采用無滑移壁面條件，計算的邊界條件針對不同的游動問題進(jìn)行不同的設(shè)置，這在下文將會詳細(xì)地給出。

1.2 動網(wǎng)格方法

考慮到魚類的游動涉及復(fù)雜的剛性或柔性變形運動，在數(shù)值計算過程中仿生魚鰭或者仿生魚表面的網(wǎng)格節(jié)點隨著計算時間的遞進(jìn)不斷地發(fā)生著變化，這將改變計算域的幾何拓?fù)浣Y(jié)構(gòu)，因此需要改變計算域中的網(wǎng)格拓?fù)浣Y(jié)構(gòu)以滿足1.1 節(jié)的CFD 計算要求。在本文的研究中采用基于彈性光順的網(wǎng)格變形以及基于尺度和歪斜率標(biāo)準(zhǔn)的局部網(wǎng)格重構(gòu)的動網(wǎng)格方法。

在基于彈性光順的網(wǎng)格變形方法中，將任意兩個相鄰的網(wǎng)格節(jié)點之間的連接理想化為相互連接的彈簧。在給定的邊節(jié)點(計算物體網(wǎng)格節(jié)點)上的位移將產(chǎn)生一個力，這個力與沿連接到該節(jié)點的所有彈簧的位移成比例。利用胡克定律，作用于網(wǎng)格節(jié)點i上的力Fi可以寫為

式中，ni為連接于網(wǎng)格節(jié)點i的相鄰節(jié)點數(shù)量，kij為網(wǎng)格節(jié)點i與相鄰節(jié)點j之間的彈簧常數(shù)，Δxi和Δxj為網(wǎng)格節(jié)點i與相鄰節(jié)點j的位移。kij表示為

式中，kfa為彈簧常數(shù)系數(shù)。

在平衡時，連接到節(jié)點i上的所有彈簧的合力必須為零，產(chǎn)生一個迭代方程為

式中，q為迭代次數(shù)。

式（5）的收斂條件由收斂容差ε和迭代次數(shù)q控制，當(dāng)滿足ε或m時，式（5）在每個時間步長的迭代解停止。收斂容差ε為

由于移動邊界（計算物體表面）的位移已知，因此在所有內(nèi)部節(jié)點上使用雅可比掃描來求解式（5）。當(dāng)收斂條件滿足時，可得到內(nèi)部所有節(jié)點的位置：

式中，n+1和n分別表示下一個時間步和當(dāng)前時間步的位置。

當(dāng)計算物體表面的網(wǎng)格運動幅度較大時，例如胸鰭的大角度旋轉(zhuǎn)拍動運動以及身體軀干和尾鰭的大幅度柔性變形運動，基于彈性光順的網(wǎng)格變形方法不再足以保證網(wǎng)格質(zhì)量，甚至產(chǎn)生負(fù)體積。為解決這一問題，基于歪斜率和尺度標(biāo)準(zhǔn)對質(zhì)量差的網(wǎng)格單元進(jìn)行聚類和局部重構(gòu)。尺寸標(biāo)準(zhǔn)是根據(jù)最小長度尺度和最大長度尺度指定的。長度尺度小于最小長度尺度或者大于最大長度尺度以及歪斜率大于設(shè)定的歪斜率的網(wǎng)格單元格被標(biāo)記進(jìn)行重構(gòu)。整個過程如下：

（1）由歪斜率和尺度標(biāo)準(zhǔn)對網(wǎng)格單元進(jìn)行標(biāo)記；

（2）刪除被標(biāo)記的網(wǎng)格單元；

（3）在被刪除的網(wǎng)格單元的空間區(qū)域生成扭曲度優(yōu)于先前網(wǎng)格的新網(wǎng)格單元；

（4）在重構(gòu)后的網(wǎng)格單元進(jìn)行物理量的差值計算，保證與之前的物理量相同。

2 仿魚類游動數(shù)值計算方法

2.1 數(shù)值計算網(wǎng)格劃分策略

正如1.2 節(jié)所述，針對仿魚類運動體的水動力數(shù)值模擬之前，需要劃分適用于動網(wǎng)格方法的數(shù)值計算網(wǎng)格，因此應(yīng)針對魚類游動的特點，建立合理的數(shù)值計算網(wǎng)格劃分策略，使得生成的網(wǎng)格能夠滿足動網(wǎng)格方法的要求，在保證計算精度的同時減少網(wǎng)格數(shù)量。以三維矩形水翼為對象，給出適用于仿生魚水動力學(xué)數(shù)值計算的網(wǎng)格劃分策略，如圖1 所示。整個計算域為一個長方體，水翼上下、左右對稱布置于計算域內(nèi)，與計算域下游邊界的距離大于到上游邊界的距離，比值通常為3～5。將整個計算域劃分為內(nèi)外兩個域，其中內(nèi)部計算域包含著整個水翼，通常也為長方體，外部計算域為整個計算域減去內(nèi)部計算域剩余的幾何空間。內(nèi)部計算域的尺寸應(yīng)保證水翼在運動過程中不會與其邊界發(fā)生接觸。在確定了計算域的幾何拓?fù)浣Y(jié)構(gòu)之后，接下來進(jìn)行網(wǎng)格的劃分。

圖1 水翼的數(shù)值計算網(wǎng)格Fig.1 Numerical mesh of the hydrofoil

首先對內(nèi)部計算域進(jìn)行網(wǎng)格的劃分，在水翼表面劃分三角形網(wǎng)格，通過控制網(wǎng)格的稀疏程度對曲率小的區(qū)域的幾何特征進(jìn)行有效的捕捉，在內(nèi)部計算域的六個表面劃分均勻的四邊形網(wǎng)格，以水翼和內(nèi)部計算域邊界的面網(wǎng)格為基礎(chǔ)，采用Delaunay三角剖分法對內(nèi)部計算域進(jìn)行四面體非結(jié)構(gòu)網(wǎng)格劃分，其中內(nèi)部計算域邊界處的網(wǎng)格為金字塔網(wǎng)格。

之后，采用分塊結(jié)構(gòu)網(wǎng)格劃分的方法對外部計算域進(jìn)行網(wǎng)格劃分。首先，使用等比數(shù)列函數(shù)將內(nèi)部計算域六個面上均勻的四邊形網(wǎng)格拉伸到外部域的邊界生成體網(wǎng)格，控制合理的增長率以保證網(wǎng)格對流場細(xì)節(jié)的捕捉。在寬度和高度方向上網(wǎng)格迅速拉伸；在縱向上，計算物體上游區(qū)域網(wǎng)格的拉伸速度較快，而在下游的尾流區(qū)域，網(wǎng)格的增長率保持在2%以下，以保持較高的尾渦分辨率。最后，對于外部計算域的其他區(qū)域，以周圍的表面網(wǎng)格為條件直接切分生成六面體網(wǎng)格。整個計算域的網(wǎng)格劃分情況如圖1所示，在本文后續(xù)的數(shù)值計算中均采用這種網(wǎng)格劃分策略。

2.2 表面網(wǎng)格更新方法

自然界中的魚類在游動過程中，身體上的各類鰭及身體軀干隨著時間的推進(jìn)不斷地發(fā)生著變化，因此在數(shù)值計算過程中它們的幾何形狀或者空間位置同樣需要發(fā)生不斷的變化，具體表現(xiàn)為表面網(wǎng)格節(jié)點坐標(biāo)的變化。FLUENT 的動網(wǎng)格功能提供了強大的二次開發(fā)能力，研究者可以基于FLUENT的用戶自定義函數(shù)（UDF）功能，采用C 語言編寫二次開發(fā)程序，對數(shù)值計算模型的剛性或者柔性運動進(jìn)行控制。

針對動網(wǎng)格的UDF 程序的開發(fā)，F(xiàn)LUENT 提供了多種宏命令，本文的數(shù)值計算中，借助DE?FINE_CG_MOTION 宏命令實現(xiàn)各種復(fù)雜的剛性體運動，例如拍動翼、尾鰭以及胸鰭的剛性運動；借助DEFINE_GRID_MOTION 宏命令實現(xiàn)各種柔性變形運動，例如身體軀干以及各類鰭的柔性運動。采用DEFINE_CG_MOTION 對采用剛性運動的胸鰭或尾鰭等魚鰭進(jìn)行控制時，可以通過控制魚鰭的三個線速度和角速度實現(xiàn)對魚鰭表面網(wǎng)格節(jié)點的三維坐標(biāo)進(jìn)行整體更新，具體過程如圖2 所示。采用DEFINE_GRID_MOTION 對仿生魚的柔性運動進(jìn)行控制時，首先對面網(wǎng)格單元的所有節(jié)點進(jìn)行遍歷，在遍歷過程中，首先讀取當(dāng)前時間步下仿生魚表面網(wǎng)格節(jié)點的三維坐標(biāo)，再通過給定的運動模型計算得到下一時刻仿生魚表面網(wǎng)格節(jié)點的坐標(biāo)，在時間步更新之后對仿生魚表面網(wǎng)格節(jié)點的三維坐標(biāo)進(jìn)行更新，具體過程如圖3所示。

圖2 剛性運動求解過程Fig.2 Solving process of rigid motion

最后，由1.2 節(jié)的動網(wǎng)格方法對計算域網(wǎng)格進(jìn)行變形與重構(gòu)，由1.1 節(jié)的數(shù)值方法對流場進(jìn)行求解，得到魚鰭或身體軀干受到的水動力。涉及自主游動時，在得到水動力之后需要運行魚體-流體耦合求解程序，對仿生魚的游動速度進(jìn)行求解。

2.3 魚體-流體耦合求解方法

魚類依靠身體上的各類魚鰭能夠進(jìn)行各種各樣的游動，然而依靠身體和尾鰭的游動通?？梢钥醋魇瞧矫鎯?nèi)的游動，這時尾鰭和身體軀干的運動通常為平面內(nèi)的運動。在平面內(nèi)的運動包含兩個平移運動和一個旋轉(zhuǎn)運動，考慮到魚類的向前游動性能是水下航行器所追求的，因此相關(guān)學(xué)者對平面的游動進(jìn)行進(jìn)一步的簡化，只保留前進(jìn)方向上的運動并進(jìn)行研究。針對魚類的自主游動數(shù)值模擬，通常不考慮魚體柔性變形過程中的內(nèi)力，僅考慮魚體在外界水動力作用下的運動[9-13]。對于平面內(nèi)的運動，由質(zhì)心運動定理：

式中，m為仿生魚的質(zhì)量，UX和UY分別為XOY平面下的兩個線速度，aX和aY分別為XOY平面下的兩個線加速度，F(xiàn)X和FY分別為XOY平面下的兩個水動力，ωz、αz、Iz和Mz分別為質(zhì)心坐標(biāo)系下的角速度、角加速度、轉(zhuǎn)動慣量和水動力矩。FX、FY和Mz分別通過對仿生魚表面的粘性力和壓力進(jìn)行積分得到：

式中，p為局部面元dS的壓力矢量，nj為j方向單位矢量，τij為黏性應(yīng)力張量，r為相對于質(zhì)心的位置矢量。

將式（9）計算得到的水動力和力矩代入式（8），即可對仿生魚的運動進(jìn)行求解。針對這一問題的求解通常分為兩種方法：一種是求解式（8）得到仿生魚的線加速度和角加速度，再對加速度進(jìn)行積分得到游動速度；另一種是由式（8）進(jìn)行直接求解得到仿生魚的游動速度。這兩種方法在數(shù)值求解精度上均能夠滿足數(shù)值計算的需要，在下文的數(shù)值計算中將對這兩種方法計算得到的游動速度進(jìn)行對比。兩種方法的計算流程如下：

方法1：

（1）初始時刻，仿生魚的速度和加速度均為零；

（2）第n個時間步，將水動力和力矩代入式（8），計算得到仿生魚的線加速度和角加速度；

（3）采用二階龍格庫塔積分獲得第n個時間步迭代結(jié)束時仿生魚的線速度和角速度，

（4）時間步更新，進(jìn)入下一次求解。

方法2：

（1）初始時刻，仿生魚的速度和加速度均為零；

（2）第n個時間步，將水動力和力矩代入式（8），采用三階向后差分方法，直接計算獲得第n個時間步迭代結(jié)束之后仿生魚的線速度和角速度，

（3）時間步更新，進(jìn)入下一次求解。

在得到仿生魚的速度信息之后，通過無滑移壁面邊界條件將仿生魚表面的速度施加于流場。無滑移壁面邊界條件為

式中，U為仿生魚表面流體的速度，Us為仿生魚表面的速度。

式中：UT=(UX,UY)為仿生魚的平移速度；UR=ωz×r為仿生魚的旋轉(zhuǎn)角速度；UDEF為仿生魚表面給定的柔性變形速度，由給定的動力學(xué)模型得到。隨后進(jìn)行下一時間步的數(shù)值計算，通過時間步的不斷遞進(jìn)，實現(xiàn)仿生魚在自身推進(jìn)下的自主游動。

3 數(shù)值計算案例

3.1 單個鰭在均勻來流中的水動力計算

魚類身體上有著多個種類的鰭，每種鰭都有著特定的功能，例如尾鰭主要進(jìn)行快速推進(jìn)，胸鰭主要進(jìn)行靈活的操縱。因此，為了便于計算以及機理分析，學(xué)者們以單個鰭為研究對象，對水動力性能及機理進(jìn)行數(shù)值研究，以求了解魚類的游動機理。尾鰭的高效推進(jìn)性能和胸鰭的靈活操縱性能是研究者所密切關(guān)心的，因而得到較為廣泛的研究，并形成了若干經(jīng)典的計算模型。

（1）仿尾鰭運動模型

首先以尾鰭為例，學(xué)者們通常以作升沉-俯仰運動的拍動翼模擬尾鰭的運動，并且已經(jīng)成為尾鰭的水動力數(shù)值研究的經(jīng)典驗證算例。Dong[14]以三維矩形翼為研究對象，針對水翼在不同運動參數(shù)下的推進(jìn)性能開展了詳細(xì)的水池實驗，水翼在均勻來流Uin中升沉-俯仰運動，水翼的展長Sfin和弦長Cfin分別為0.6 m 和0.1 m，水翼數(shù)值計算模型如圖4 所示。關(guān)于水翼的詳細(xì)實驗研究見文獻(xiàn)[14]，本文不做過多描述，只對水翼的運動模型和相關(guān)運動參數(shù)進(jìn)行介紹。水翼的升沉運動hfin(t)和俯仰運動θfin(t)為

圖4 水翼數(shù)值計算模型Fig.4 Computational model of the hydrofoil

式中，h0和θ0分別為水翼的升沉運動幅值和俯仰運動幅值，f為運動頻率，?為俯仰運動與升沉運動的相位差。水翼在運動時的瞬時有效攻角α(t)定義為

水翼在作升沉-俯仰耦合運動時，影響水翼推進(jìn)性能的一個關(guān)鍵無因次參數(shù)為斯特勞哈爾數(shù)St，定義為

水翼在給定的運動參數(shù)下進(jìn)行運動時的推力系數(shù)CT定義為

式中，Tfin為水翼一個運動周內(nèi)推力的均值。

在使用CFD 方法對流體中的物體進(jìn)行水動力數(shù)值計算時，用于數(shù)值計算的網(wǎng)格的合適與否直接影響著數(shù)值計算結(jié)果的精度。采用前文給出的數(shù)值計算方法和網(wǎng)格劃分策略，其中數(shù)值計算網(wǎng)格與圖1 給出的一致，計算域的出口設(shè)置為壓力出口，計算域的其他面均為速度入口。由2.1 節(jié)給出的網(wǎng)格劃分策略，對內(nèi)部計算域和外部計算域劃分三種不同數(shù)量的網(wǎng)格，按照網(wǎng)格由密集到稀疏的順序分別命名為mesh 1、mesh 2 和mesh 3。采用這三種網(wǎng)格對水翼在運動參數(shù)為h0/Cfin=0.8、St=0.5、?=90°、Uin=0.4 m/s和θ0=20°下的推進(jìn)性能進(jìn)行數(shù)值模擬，由計算得到的水翼的推力系數(shù)CT對網(wǎng)格數(shù)量的無關(guān)性進(jìn)行分析，如表1所示。

表1 不同網(wǎng)格計算值Tab.1 Computation by different meshes

基于Stern 等[15]提出基于Richard 外推（Richard extrapolation,RE）法給出網(wǎng)格數(shù)量的收斂性驗證的方法和步驟。在進(jìn)行網(wǎng)格無關(guān)性驗證時，首先需要給出網(wǎng)格細(xì)化率rg。一般來說，rg的定義是粗細(xì)網(wǎng)格的單元格大小之比。Baek 等[16]給出了采用非結(jié)構(gòu)網(wǎng)格進(jìn)行CFD 數(shù)值模擬時網(wǎng)格細(xì)化率的計算方法，這個方法已經(jīng)被廣泛地應(yīng)用于船舶螺旋槳、水下航行器以及水面船等多種涉及CFD網(wǎng)格數(shù)量的無關(guān)性驗證中[17]。rg的計算公式為

式中，Nfine和Ncoarse分別表示相鄰尺度的網(wǎng)格系統(tǒng)中細(xì)密和粗糙的網(wǎng)格的總數(shù)目，d表示求解問題的空間維度。rg在本節(jié)的驗證中約為1.2。

以si(i=1,2,3)分別表示mesh 1、mesh 2 和mesh 3 計算得到的結(jié)果，則不同網(wǎng)格系統(tǒng)計算結(jié)果的差值為

網(wǎng)格系統(tǒng)的收斂率Rg定義為

Stern等[15]依據(jù)Rg的值將網(wǎng)格系統(tǒng)的收斂情況分為以下四種：

（a）0＜Rg＜1時，單調(diào)收斂；

（b）Rg＜0且 ||

Rg＜1時，振蕩收斂；

（c）Rg＞1時，單調(diào)發(fā)散；

（d）Rg＜0且 ||

Rg＞1時，振蕩發(fā)散。

三組網(wǎng)格的Rg約為0.104，表明本文所采用的網(wǎng)格劃分策略得到的網(wǎng)格在求解仿魚類運動體的水動力性能時具有較好的收斂性。在本文后續(xù)的研究中所采用的數(shù)值計算網(wǎng)格均采用這種方法進(jìn)行網(wǎng)格的獨立驗證。

水翼的運動參數(shù)為h0/Cfin=1.0，St以0.05為間隔從0.25變化至0.4，?=90°，Uin=0.4 m/s，α(t)的最大值等于30°，獲得每個St下水翼的CT。圖5顯示了本文計算結(jié)果與Dong[14]的實驗結(jié)果和數(shù)值計算結(jié)果的對比，從圖5 中可以看出本文的計算結(jié)果與實驗結(jié)果吻合較好，并且優(yōu)于Dong[14]采用的前緣渦分離（LVS）的BEM 所計算得到的結(jié)果。由于BEM 不考慮流體的粘性，無法對水翼受到流體作用的粘性阻力進(jìn)行精確的計算，因而得到的推力系數(shù)較實驗結(jié)果和本文的計算結(jié)果偏高。從圖中可以看出，本文基于CFD技術(shù)建立的仿生水動力學(xué)求解方法能夠?qū)Ψ律膭右淼乃畡恿π阅苓M(jìn)行精確的求解。本文采用Hunt 等[18]提出的Q準(zhǔn)則對流場中的三維渦結(jié)構(gòu)進(jìn)行可視化處理，通過Q等值面對流場中的三維渦結(jié)構(gòu)進(jìn)行顯示，Q的定義為Q=(‖W‖2-‖S‖2)/2，式中W和S分別表示速度梯度的反對稱張量和對稱張量，‖ ‖為矩陣的歐幾里德范數(shù)。根據(jù)Hunt等[18]的定義，當(dāng)Q大于零時，Q等值面為旋渦結(jié)構(gòu)，本文中對流場的三維渦結(jié)構(gòu)的顯示均采用Q準(zhǔn)則。圖6 給出了水翼在運動過程中某一瞬間流場的三維渦結(jié)構(gòu)。從圖中可以看出水翼上表面靠近前緣的區(qū)域產(chǎn)生了一個強度較大的前緣渦，水翼上表面被前緣渦附著的區(qū)域被強度較大的低壓所覆蓋，如圖7 所示；在水翼的兩個梢部有明顯的梢渦產(chǎn)生，兩個梢渦向上卷曲，與前緣渦之間有著明顯的分界；在水翼的尾緣處存在明顯的尾緣渦，尾緣渦與水翼的尾緣相連并脫落于尾流場中。由三維渦結(jié)構(gòu)特征可以看出，本文采用非結(jié)構(gòu)網(wǎng)格與結(jié)構(gòu)網(wǎng)格混合的分塊網(wǎng)格劃分方法，對包含計算物體的流場區(qū)域的加密能夠有效地捕捉水翼表面產(chǎn)生并脫落于流場中的渦結(jié)構(gòu)特征。

圖6 水翼產(chǎn)生的三維渦結(jié)構(gòu)Fig.6 Three-dimensional vortex structure generated by hydrofoil

圖7 水翼表面的壓力分布云圖Fig.7 Pressure distribution on hydrofoil surface

（2）仿胸鰭運動模型

接下來以胸鰭為例，魚類的胸鰭有著多個旋轉(zhuǎn)運動組合而成的復(fù)雜運動，與尾鰭的運動模式有著明顯的區(qū)別，基于前人對魚類胸鰭的實驗研究中給出的幾何模型作為數(shù)值計算模型，對本文基于CFD技術(shù)建立的仿生水動力學(xué)求解方法在求解胸鰭水動力性能上的有效性進(jìn)行驗證。Kato和Furushima[19]針對黑鱸魚的實驗研究建立了真實黑鱸魚胸鰭的幾何模型，如圖8所示。與Suzuki等[20]對黑鱸魚胸鰭的水動力實驗結(jié)果進(jìn)行對比，在Suzuki 等[20]的實驗中，胸鰭的展長Spec為0.1 m，弦長Cpec為0.08 m，厚度為0.7 mm，展弦比（AR=，Apec=0.005 834 m2為胸鰭的投影面積）約為1.7?？拷~體的一側(cè)為鰭根，遠(yuǎn)離魚體的一側(cè)為鰭梢，迎流的一側(cè)為前緣，去流的一側(cè)為尾緣，鰭根與轉(zhuǎn)動中心的距離為0.055 m。

圖8 胸鰭數(shù)值計算模型Fig.8 Computational model of the pectoral

胸鰭的運動由前后的拍動運動和縱傾運動組成，胸鰭的拍動運動和縱傾運動如圖9 所示。胸鰭的拍動運動和縱傾運動方程[20]為

圖9 胸鰭的旋轉(zhuǎn)運動Fig.9 Rotational motion of the pectoral fin

式中，θRC為拍動角的周期平均值，θRA為拍動角幅值，ωpec=2πf為胸鰭運動的圓頻率，θFC為縱傾角的周期平均值，θFA為縱傾角幅值，ΔθF為縱傾運動與拍動運動的相位差。

胸鰭在均勻來流Uin=0.201 m/s 中進(jìn)行運動，運動參數(shù)設(shè)置為：θRC=30°，θRA=30°，θFC=30°，θFA=30°，ΔθF=90°，f=2.0 Hz。通過數(shù)值計算獲得胸鰭的縱向力系數(shù)CXpec和側(cè)向力系數(shù)CYpec，與Suzuki 等[20]的實驗結(jié)果進(jìn)行對比，驗證本文的數(shù)值計算方法在求解胸鰭水動力性能上的有效性，縱向力系數(shù)CXpec和側(cè)向力系數(shù)CYpec定義為

式中，F(xiàn)Xpec為胸鰭的縱向力，F(xiàn)Ypec為胸鰭的側(cè)向力。

采用前文給出的網(wǎng)格劃分策略生成胸鰭的數(shù)值計算網(wǎng)格，如圖10 所示，整個計算域的網(wǎng)格數(shù)量為3 335 564，采用與水翼相同的數(shù)值計算方法。

圖11給出了本文的數(shù)值計算方法得到的胸鰭縱向力系數(shù)CXpec和側(cè)向力系數(shù)CYpec與Suzuki 等[20]的實驗結(jié)果的對比。從圖中可以看出，實驗值的變化趨勢發(fā)生改變時會出現(xiàn)小幅的波動，并且計算結(jié)果與實驗值的幅值存在一定的差異，但是兩者的總體變化趨勢是一致的，并且吻合得較好。由以上分析可以得出，本文基于CFD 技術(shù)建立的仿生水動力學(xué)求解方法能夠?qū)哂卸嘈D(zhuǎn)運動的胸鰭的水動力性能進(jìn)行有效的求解。

圖11 胸鰭的水動力系數(shù)對比Fig.11 Comparison of hydrodynamic coefficients of pectoral fin

圖12 和圖13 分別給出了胸鰭運動過程t/T=0.875 時流場的三維渦結(jié)構(gòu)和表面壓力分布云圖。從圖12 可以看出在胸鰭的前緣附著一個強度較大的前緣渦，這導(dǎo)致了一個強度較大的低壓區(qū)域（圖13）；在胸鰭的梢部、根部以及尾緣處均產(chǎn)生了強度較大的渦，并且這些渦與胸鰭連接在一起并脫落于尾流場中，在尾流場中存在著先前運動周期所產(chǎn)生的渦環(huán)以及梢渦。盡管胸鰭的運動為空間內(nèi)復(fù)雜的旋轉(zhuǎn)運動，會給網(wǎng)格帶來高度的扭曲變形，但是本文的網(wǎng)格劃分策略同樣能夠?qū)α鲌鲋袦u結(jié)構(gòu)進(jìn)行較好的捕捉。胸鰭上表面的前緣區(qū)域附著一個強度較大的前緣渦，前緣渦附著的區(qū)域被強度較大的低壓所覆蓋。

圖12 胸鰭產(chǎn)生的三維渦結(jié)構(gòu)Fig.12 Three-dimensional vortex structure generated by pectoral fin

圖13 胸鰭表面的壓力分布云圖Fig.13 Pressure distribution on pectoral fin surface

3.2 仿生魚在均勻來流中的水動力計算

魚類在游動過程中身體上的各類鰭必定與軀體產(chǎn)生相互干擾，采用單個鰭進(jìn)行魚類游動機理的研究無法計入這種干擾的影響。因此以整魚為數(shù)值計算模型，研究各類鰭在軀體干擾下的水動力性能對進(jìn)一步掌握魚類游動機理是十分有必要的。整魚的運動涉及軀體的柔性波動運動，因此難以開展精確的實驗研究，目前針對這類問題的數(shù)值計算方法驗證通常以3.1 節(jié)中的水翼或者胸鰭作為替代，這不能充分地證明計算方法是合理有效的。為有效開展仿生魚在均勻來流中的水動力性能數(shù)值驗證，本節(jié)以Zhu 等[7]給出的RoboTuna 仿生機器魚幾何模型為驗證算例。國內(nèi)外多個科研團(tuán)隊以Ro?boTuna為數(shù)值計算模型，采用BEM 或者CFD 方法對其游動過程中的水動力性能進(jìn)行了數(shù)值計算及試驗驗證，通過與這些數(shù)值計算結(jié)果的對比，能夠有效地驗證本文所采用的數(shù)值計算方法在求解整魚水動力性能上的有效性。

仿生魚數(shù)值計算模型如圖14 所示，圖中紅色虛線表示魚體的中線，O-XYZ為全局坐標(biāo)系，o-xyz為固定在魚體上的隨體坐標(biāo)系，x軸與魚體的中線重合由頭部指向尾鰭，坐標(biāo)原點位于距離魚體前端30%身體軀干（不包含尾鰭）長度處，y軸指向魚體右側(cè)，z軸滿足右手準(zhǔn)則。仿生魚幾何外形上下左右對稱，身體軀干的長度L為1 m，外形為紡錘形，垂直于x軸的截面為橢圓形，橢圓的長短軸之比為3:2，長軸與z軸平行。因此，由身體軀干在中縱剖面上的輪廓線可以確定身體軀干的幾何外形，輪廓線的曲線方程為

圖14 仿生魚幾何模型Fig.14 Geometric model of the bionic fish

尾鰭前緣x(z)LE和尾緣x(z)TE的曲線方程分別為

仿生魚在游動過程中，身體軀干進(jìn)行側(cè)向的柔性波動運動，尾鰭跟隨身體軀干進(jìn)行側(cè)向平移運動的同時圍繞著尾柄末端進(jìn)行擺動運動。身體軀干的側(cè)向柔性波動運動方程為

式中：a(x)為柔性波動運動的振幅包絡(luò)線函數(shù)；k為波數(shù)，定義為k=2π/λ，λ為波長；a1和a2均為常數(shù)，取不同的值可以調(diào)整尾鰭的側(cè)向平移幅值。

尾鰭在尾柄末端的帶動下作側(cè)向平移運動的同時繞著尾柄末端作擺動運動，尾鰭的擺動運動方程為

式中，xp為尾柄末端坐標(biāo)，θm為擺角幅值，Δθ為尾鰭的擺動運動與平移運動之間的相位差。

在驗證計算中，仿生魚被置于均勻的來流Uin=0.7 m/s 中，身體軀干的運動參數(shù)為a1=0.002 36，a2=-0.114，f=2.07 Hz，λ=1.675 m，尾鰭的運動參數(shù)為θm=18.76°，Δθ=85°。整個計算域網(wǎng)格數(shù)量為3 693 734，如圖15 所示。圖16 給出了仿生魚縱向力FT的時歷曲線，并與Zhu 等[7]的BEM 計算結(jié)果以及Zhang等[21]和Chang 等[22]的CFD 計算結(jié)果進(jìn)行了對比。從圖16 中可以看出，本文的計算結(jié)果與文獻(xiàn)中的計算結(jié)果吻合較好，F(xiàn)T的幅值小于BEM 的計算結(jié)果，這是因為BEM 沒有考慮流體的粘性而導(dǎo)致的。此外，本文的數(shù)值計算方法得到的結(jié)果與Chang 等[22]的計算結(jié)果有著較高的吻合度，表明本文的數(shù)值計算方法可以有效地獲得仿生魚做柔性運動時的水動力。

圖15 數(shù)值計算網(wǎng)格Fig.15 Numerical mesh

圖16 仿生魚推力曲線Fig.16 Thrust curve of the bionic fish

3.3 仿生魚在靜水中的3-DoF自主游動數(shù)值模擬

（1）鰻魚數(shù)值計算驗證

Kern 和Koumoutsakos[23]對鰻鱺游動模式的仿生魚在水平面內(nèi)的3-DoF 自主游動進(jìn)行了數(shù)值計算，并且給出了仿生魚數(shù)值計算模型的詳細(xì)幾何參數(shù)以及運動模型，針對仿生魚在給定的一組運動參數(shù)下的自主游動進(jìn)行了數(shù)值計算，并給出了最終的巡游速度，因此非常適用于仿生魚自主游動的數(shù)值驗證。仿生魚體長La=1 m，仿生魚的三維幾何模型由空間上變化的橢圓橫截面構(gòu)成，橢圓的長短軸由仿生魚高度和寬度方向上的輪廓線h(s)和w(s)控制，w(s)的定義為

式中，s為仿生魚表面到頭部頂點的縱向距離，wh=sb=0.04La，st=0.95La，wt=0.01La。h(s)的定義為

式中，a=0.51La，b=0.08La。

仿生魚寬度和高度方向上的輪廓如圖17所示。

圖17 仿生魚寬度和高度方向上的輪廓圖Fig.17 Outline of the bionic fish in the width and height directions

仿生魚在自主游動過程中某一瞬間身體的柔性波動如圖18所示。在自主游動過程中，仿生魚身體的側(cè)向柔性波動方程為

圖18 仿生魚身體幾何構(gòu)型Fig.18 Geometry of the anguilliform fish

在驗證計算中，流體的動粘滯系數(shù)μ=1.4×10-4kg/ms，流體的密度ρ和仿生魚的密度ρfish相同且均為1 kg/m3，f=1 Hz。網(wǎng)格劃分策略及數(shù)值計算方法與3.2節(jié)的一致，由于涉及仿生魚的自主游動，將入口邊界的流體速度和壓力梯度設(shè)置為零，出口邊界的速度梯度和壓力梯度設(shè)置為零。仿生魚由靜止開始擺動身體進(jìn)行游動，圖19 給出了鰻鱺游動模式下仿生魚的前進(jìn)速度Uf的時歷曲線，其中紅色實線是三次向后差分方法的計算結(jié)果，黑色點虛線是二階龍格庫塔積分的計算結(jié)果，從圖中可以看出兩種數(shù)值算法的計算結(jié)果是一致的。本文的數(shù)值計算方法得到仿生魚最終的巡游速度Ucf=0.381La/s，與Kern和Koumoutsakos等[23]的數(shù)值計算中得到的0.40La/s巡游速度吻合較好（相差-4.75%），表明本文的數(shù)值計算方法能夠?qū)Ψ律~的自主游動過程進(jìn)行精確的模擬。圖20 給出了鰻魚在自主游動過程中某一瞬間流場的三維渦結(jié)構(gòu)。從圖中可以看出鰻魚的尾渦為雙排交錯排列結(jié)構(gòu)，漩渦在身體彎曲的一側(cè)產(chǎn)生，隨著身體上行進(jìn)波向后的傳播，逐漸由頭部向尾鰭移動，在身體軀干的后半段形成一個明顯的“發(fā)卡”渦，隨著尾鰭運動方向的轉(zhuǎn)變，脫落于為流場中，形成一個一個尾渦結(jié)構(gòu)。

圖19 前進(jìn)速度Fig.19 Forward velocity

圖20 鰻魚產(chǎn)生的三維渦結(jié)構(gòu)Fig.20 Three-dimensional vortex structures generated by the eel

（2）金槍魚數(shù)值計算驗證

接下來以哈爾濱工程大學(xué)水下機器人技術(shù)重點實驗室研制的“仿生-I”號機器魚的定向游動實驗，進(jìn)一步地驗證本文的數(shù)值計算方法在求解仿生魚自主游動上的有效性?！胺律?I”號機器魚及數(shù)值計算模型如圖21 所示。關(guān)于“仿生-I”號機器魚的詳細(xì)試驗見文獻(xiàn)[24，25]，關(guān)于仿生魚的運動模型以及運動參數(shù)見文獻(xiàn)[12]，文獻(xiàn)[26]給出了仿生魚數(shù)值計算模型的詳細(xì)型值，包括身體軀干、第一背鰭、第二背鰭、胸鰭、臀鰭和尾鰭，可以用于數(shù)值計算模型的幾何建立，以上內(nèi)容本文不再給出。通過改變仿生魚的游動頻率f，得到最終的巡游速度Ucf，圖22 顯示了由本文3-DoF 數(shù)值計算方法得到的結(jié)果與“仿生-I”號機器魚的試驗結(jié)果[22]以及蘇玉民等[27]的1-DoF 數(shù)值計算結(jié)果對比，圖中虛心點為計算或試驗的結(jié)果，點畫線為二次擬合的結(jié)果，從圖中可以看出本文的數(shù)值計算結(jié)果與試驗結(jié)果的趨勢一致，Ucf隨f的增加均表現(xiàn)為非線性增長趨勢，并且優(yōu)于蘇玉民等[27]的1-DoF 數(shù)值計算結(jié)果中Ucf隨f的增加表現(xiàn)為線性增長趨勢。另外，除了f為0.8 Hz時本文的計算結(jié)果與試驗值相差較大（20.5%）以外，其他情況下的差值均在10%以下。以上分析進(jìn)一步驗證了本文的數(shù)值計算方法在求解魚類自主游動方面的有效性。

圖21“仿生-I”號機器魚及數(shù)值計算模型Fig.21“FangSheng-I”robotic fish and computational model

圖22“仿生-I”號機器魚的巡游實驗對比Fig.22 Comparison of the cruising experi?ments of“Fang Sheng-I”robotic fish

4 結(jié) 論

仿生水動力學(xué)長期都是計算流體力學(xué)中的研究熱點，其關(guān)注的是水中生物高效、快速和高機動游動行為的機理，在水中推進(jìn)和減阻方面有著重要的工程應(yīng)用前景。本文基于CFD方法，結(jié)合動網(wǎng)格方法，給出了仿生水動力學(xué)數(shù)值計算方法，針對單個魚鰭的水動力、均勻來流中仿生魚的水動力以及靜水中仿生魚的自主游動等仿生水動力學(xué)問題開展了研究，并與文獻(xiàn)中經(jīng)典的且可對比的實驗或數(shù)值結(jié)果進(jìn)行了比較，驗證了本文的數(shù)值計算方法在求解仿生水動力學(xué)問題上的通用性和有效性。本文的主要結(jié)論如下：

（1）基于彈性光順和局部網(wǎng)格重構(gòu)的動網(wǎng)格方法能夠處理仿生尾鰭的大幅拍動運動、胸鰭的大角度旋轉(zhuǎn)運動以及身體軀干的柔性波動運動，在處理各類仿生運動帶來的網(wǎng)格大變形運動上有著較強的健壯性；

（2）采用多塊網(wǎng)格劃分策略能夠有效地處理仿生水動力學(xué)的計算問題，對包含運動體的內(nèi)部區(qū)域采用非結(jié)構(gòu)網(wǎng)格劃分方法，不包含運動體的外部區(qū)域采用結(jié)構(gòu)網(wǎng)格劃分方法，這種網(wǎng)格劃分策略能夠在保證計算精度的前提下，減少網(wǎng)格數(shù)量，提高計算效率；拍動翼、胸鰭以及整魚的尾流場中的三維渦結(jié)構(gòu)的顯示結(jié)果表明，這種網(wǎng)格劃分方法能夠?qū)α鲌鲋泻诵牡臏u結(jié)構(gòu)特征進(jìn)行捕捉。

（3）通過與文獻(xiàn)中經(jīng)典的拍動翼和胸鰭的水動力試驗結(jié)果以及仿生魚在均勻來流中進(jìn)行約束游動時的水動力性能的對比，以及與文獻(xiàn)中的仿生魚在均勻來流中的約束游動或者靜水中的自主游動數(shù)值計算的對比，表明基于CFD方法的二次開發(fā)能夠?qū)Ψ律畡恿W(xué)問題進(jìn)行有效的求解。

本文的數(shù)值計算方法為仿魚類游動水動力學(xué)走上工程應(yīng)用提供了有力的手段，能夠服務(wù)于仿生推進(jìn)器以及魚類游動的機理研究，本領(lǐng)域的研究者同樣能夠依據(jù)本文給出的各種經(jīng)典驗證算例對不同的仿生水動力學(xué)數(shù)值計算方法開展驗證。