基于CFD的仿生水動力學數值計算方法及其驗證

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

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

0 引 言

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

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

1 CFD數值計算方法

1.1 控制方程

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

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

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

1.2 動網格方法

考慮到魚類的游動涉及復雜的剛性或柔性變形運動，在數值計算過程中仿生魚鰭或者仿生魚表面的網格節點隨著計算時間的遞進不斷地發生著變化，這將改變計算域的幾何拓撲結構，因此需要改變計算域中的網格拓撲結構以滿足1.1 節的CFD 計算要求。在本文的研究中采用基于彈性光順的網格變形以及基于尺度和歪斜率標準的局部網格重構的動網格方法。

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

式中，ni為連接于網格節點i的相鄰節點數量，kij為網格節點i與相鄰節點j之間的彈簧常數，Δxi和Δxj為網格節點i與相鄰節點j的位移。kij表示為

式中，kfa為彈簧常數系數。

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

式中，q為迭代次數。

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

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

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

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

（1）由歪斜率和尺度標準對網格單元進行標記；

（2）刪除被標記的網格單元；

（3）在被刪除的網格單元的空間區域生成扭曲度優于先前網格的新網格單元；

（4）在重構后的網格單元進行物理量的差值計算，保證與之前的物理量相同。

2 仿魚類游動數值計算方法

2.1 數值計算網格劃分策略

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

圖1 水翼的數值計算網格Fig.1 Numerical mesh of the hydrofoil

首先對內部計算域進行網格的劃分，在水翼表面劃分三角形網格，通過控制網格的稀疏程度對曲率小的區域的幾何特征進行有效的捕捉，在內部計算域的六個表面劃分均勻的四邊形網格，以水翼和內部計算域邊界的面網格為基礎，采用Delaunay三角剖分法對內部計算域進行四面體非結構網格劃分，其中內部計算域邊界處的網格為金字塔網格。

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

2.2 表面網格更新方法

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

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

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

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

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

魚類依靠身體上的各類魚鰭能夠進行各種各樣的游動，然而依靠身體和尾鰭的游動通常可以看作是平面內的游動，這時尾鰭和身體軀干的運動通常為平面內的運動。在平面內的運動包含兩個平移運動和一個旋轉運動，考慮到魚類的向前游動性能是水下航行器所追求的，因此相關學者對平面的游動進行進一步的簡化，只保留前進方向上的運動并進行研究。針對魚類的自主游動數值模擬，通常不考慮魚體柔性變形過程中的內力，僅考慮魚體在外界水動力作用下的運動[9-13]。對于平面內的運動，由質心運動定理：

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

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

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

方法1：

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

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

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

（4）時間步更新，進入下一次求解。

方法2：

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

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

（3）時間步更新，進入下一次求解。

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

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

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

3 數值計算案例

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

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

（1）仿尾鰭運動模型

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

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

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

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

水翼在給定的運動參數下進行運動時的推力系數CT定義為

式中，Tfin為水翼一個運動周內推力的均值。

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

表1 不同網格計算值Tab.1 Computation by different meshes

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

式中，Nfine和Ncoarse分別表示相鄰尺度的網格系統中細密和粗糙的網格的總數目，d表示求解問題的空間維度。rg在本節的驗證中約為1.2。

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

網格系統的收斂率Rg定義為

Stern等[15]依據Rg的值將網格系統的收斂情況分為以下四種：

（a）0＜Rg＜1時，單調收斂；

（b）Rg＜0且 ||

Rg＜1時，振蕩收斂；

（c）Rg＞1時，單調發散；

（d）Rg＜0且 ||

Rg＞1時，振蕩發散。

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

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

圖6 水翼產生的三維渦結構Fig.6 Three-dimensional vortex structure generated by hydrofoil

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

（2）仿胸鰭運動模型

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

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

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

圖9 胸鰭的旋轉運動Fig.9 Rotational motion of the pectoral fin

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

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

式中，FXpec為胸鰭的縱向力，FYpec為胸鰭的側向力。

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

圖11給出了本文的數值計算方法得到的胸鰭縱向力系數CXpec和側向力系數CYpec與Suzuki 等[20]的實驗結果的對比。從圖中可以看出，實驗值的變化趨勢發生改變時會出現小幅的波動，并且計算結果與實驗值的幅值存在一定的差異，但是兩者的總體變化趨勢是一致的，并且吻合得較好。由以上分析可以得出，本文基于CFD 技術建立的仿生水動力學求解方法能夠對具有多旋轉運動的胸鰭的水動力性能進行有效的求解。

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

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

圖12 胸鰭產生的三維渦結構Fig.12 Three-dimensional vortex structure generated by pectoral fin

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

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

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

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

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

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

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

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

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

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

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

圖15 數值計算網格Fig.15 Numerical mesh

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

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

（1）鰻魚數值計算驗證

Kern 和Koumoutsakos[23]對鰻鱺游動模式的仿生魚在水平面內的3-DoF 自主游動進行了數值計算，并且給出了仿生魚數值計算模型的詳細幾何參數以及運動模型，針對仿生魚在給定的一組運動參數下的自主游動進行了數值計算，并給出了最終的巡游速度，因此非常適用于仿生魚自主游動的數值驗證。仿生魚體長La=1 m，仿生魚的三維幾何模型由空間上變化的橢圓橫截面構成，橢圓的長短軸由仿生魚高度和寬度方向上的輪廓線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所示。在自主游動過程中，仿生魚身體的側向柔性波動方程為

圖18 仿生魚身體幾何構型Fig.18 Geometry of the anguilliform fish

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

圖19 前進速度Fig.19 Forward velocity

圖20 鰻魚產生的三維渦結構Fig.20 Three-dimensional vortex structures generated by the eel

（2）金槍魚數值計算驗證

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

圖21“仿生-I”號機器魚及數值計算模型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 結 論

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

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

（2）采用多塊網格劃分策略能夠有效地處理仿生水動力學的計算問題，對包含運動體的內部區域采用非結構網格劃分方法，不包含運動體的外部區域采用結構網格劃分方法，這種網格劃分策略能夠在保證計算精度的前提下，減少網格數量，提高計算效率；拍動翼、胸鰭以及整魚的尾流場中的三維渦結構的顯示結果表明，這種網格劃分方法能夠對流場中核心的渦結構特征進行捕捉。

（3）通過與文獻中經典的拍動翼和胸鰭的水動力試驗結果以及仿生魚在均勻來流中進行約束游動時的水動力性能的對比，以及與文獻中的仿生魚在均勻來流中的約束游動或者靜水中的自主游動數值計算的對比，表明基于CFD方法的二次開發能夠對仿生水動力學問題進行有效的求解。

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