機器魚多鰭協同推進加速
——滑行行為數值模擬

李浩宇,李宗剛,杜亞江,王小波

(1.蘭州交通大學 機電工程學院, 蘭州 730070; 2.蘭州交通大學 機器人研究所， 蘭州 730070)

1 引言

仿生水下機器人是一種運用真實魚類的游動機理在水下進行快速游動、快速轉彎等平面或空間運動的機器人系統，因其機動性好、推進性能強以及隱身性好的優勢在水下勘探、生物監測、軍事偵察與打擊等方面具有重要的應用價值，是近年來國內外研究的熱點。

海洋中魚類的游動可分為身體和/或尾鰭(body and/or caudal fin,BCF)游動模式以及中央鰭和/或對鰭(median and/or paired fin,MPF)游動模式[1-3]，國內外科研人員通過深入研究魚類的游動機理，研發了各類獨具特色且擁有不同游動性能的仿生機器魚[4]。為分析解決機器魚運動與周圍流場結構的相互作用問題，現階段大多數研究利用數值模擬方法對流場進行處理以得到渦結構、速度矢量等物理量，進而闡明機器魚在流場中的游動機理，優化機器魚的行為。針對流場中胸鰭或尾鰭對流場的作用問題，李寧宇[5]將單胸鰭放入流場，分析研究胸鰭擺動時不同參數對流場結構的影響，發現單個胸鰭擺動時胸鰭的推力和效率達到最大時相位差應保持在90°；Zhang等[6]、Krishnadas等[7]等通過數值方法模擬了尾鰭形狀對推進性能的影響，結果表明，在相同運動規律下，仿金槍魚形狀的尾鰭在擺動時的推進性能和效率最佳；馮億坤[8]研究了機器魚在C型啟動過程中的水動力性能，分析了c型啟動2個階段產生的渦環和射流；Xu等[9]運用重疊網格法，將剛性胸/尾鰭放入流場中，僅靠胸鰭的擺動來推動魚體直游和轉彎；劉煥興等[10]、Wu等[11]對仿生機器魚的尾鰭柔性變形自主游動進行了研究。以上研究在分析單個胸鰭擺動或尾鰭波動的水動力性能中取得了一些成果，但并未考慮胸尾鰭同時運動時的協同作用。

胸尾鰭協同游動即在上述文獻的基礎上，將胸鰭的擺動和尾鰭的波動相結合，通過兩者的共同作用使魚體做出相應動作，它相較于單純的胸鰭或尾鰭推動，在推進性能與效率上有明顯優勢，通過胸鰭渦與尾鰭渦的融合，機器魚擁有更好的推進性能與機動性能；在胸尾鰭協同游動的研究中，Li等[12]分析了河豚的柔性多鰭動力學，將含柔性背鰭、臀鰭、尾鰭的仿河豚模型放入流場進行推進，分析了其水動力學性能和渦流結構。在多體動力學研究方面，Kanso等[13]分析了離散的鰻魚模型在理想流體中的自推進游動；John等[14]研究了連續的鰻魚模型在自主推進時的游動速度與尾渦的變化。

在其他仿生機器魚相關的研究工作中，謝鷗等[15]針對針對仿生機器魚作業時的環境識別問題，提出了一種基于人工側線(ALL)的近壁面流場識別方法，采用CFD方法提取數據并運用多層前饋神經網絡建立預測模型，通過驗證，該方法具有良好的預測效果；Zheng等[16]建立了由壓力傳感器陣列組成的人工側線系統來感受機器魚游動時周圍流場的變化，該研究驗證了人工側線系統的有效性和實用性；劉科顯等[17]通過理論分析、CFD仿真以及估計的方法研究了二維多關節機器魚遠場流速與攻角預測；和巖輝等[18]基于中樞模式發生器(CPG)理論結合模糊控制器提出了一種可實現機器魚定向游動的精確控制方法，通過仿真與實體實驗驗證了該算法的穩定性；梁旭[19]研制了一種高頻驅動在線便變剛度柔性仿生機器魚，并進行了該機器魚的游動特性實驗，驗證了魚類依靠改變剛度從而改變固有頻率來匹配擺動頻率的可行性；Yu等[20]研制了一種無線電控制的多關節機器魚，且通過模糊強化學習算法實現了多個機器魚之間的協調，并成功應用到國際水中機器人比賽中；此外，Yu等[21]還研究了多連桿機器魚的轉向控制問題，開發了一種使用方法控制機器魚的轉彎步態，并通過仿真和實體實驗驗證了該方法的適用性。

基于流體力學(CFD)方法，數值模擬了仿生機器魚在水平面內的胸尾鰭協同自主推進直線游動，其中胸鰭采用3自由度剛性擺動，尾鰭采用柔性波動；胸鰭的擺動采用重新設計的3自由度擺動曲線，并設置了4組不同占空比，分析占空比對機器魚游動速度以及水動力性能的影響。魚體由靜止開始，施加特定的運動規律進行直游，在經過一定時間的加速后進入穩定巡游階段。對整個過程中魚體與胸鰭的速度矢量、壓力云圖以及流場的三維渦結構進行分析，揭示魚類直線游動的機理，為之后的胸尾鰭協同推進研究提供參考。

2 機器魚的計算模型

2.1 胸鰭的計算模型

仿生機器魚的三維幾何模型與相關尺寸如圖1所示，其中胸鰭與魚體采用分離式結構，胸鰭進行3自由度擺動，兩側胸鰭初始位置貼近魚體以模仿真實魚類。該模型的運動部分有3個，分別是兩側的剛性胸鰭以及柔性尾鰭，2個胸鰭通過將上下拍動、前后拍動以及搖翼運動3個自由度運動耦合，以實現期望的運動軌跡，其中機器魚整體長度為1 000 mm，寬度為205 mm，尾鰭高度為370 mm，胸鰭展長為170 mm，最大弦長為130 mm。魚體分為PA和PB段，其中PA段為魚頭剛性部分，長度為400 mm，PB段為尾鰭柔性部分，長度為600 mm。胸鰭的運動學模型可由下式給出：

(1)

式中：φR、φFL、φF分別代表胸鰭的前后拍、上下拍、搖翼運動的歐拉角;φRC、φFC、φFLC分別代表前后拍、上下拍、搖翼運動角的平均值;φR0、φF0、φFL0分別代表3個運動的幅值; dφF、dφFL代表相位差;ω為運動角速度;t為時間。

圖1 仿生機器魚的三維模型及相關尺寸

圖2為胸鰭側視圖中隨體坐標系坐標軸與胸鰭外緣的交點隨時間變化的運動軌跡，其中每個橢圓表示特定幅值下胸鰭前后拍翼角與上下拍翼角的相對關系，所有橢圓的交點即胸鰭在t=0時刻的位置，該曲線為胸鰭前后拍與上下拍2個自由度的約束條件，此外，由于該圖為胸鰭XOZ平面的運動曲線，搖翼運動與另外2自由度的運動相對關系并不包含在該胸鰭擺動曲線內。

圖2 胸鰭的擺動曲線

“加速-滑行”游動模式是魚類等水生生物常見的游動方式，自然界中的魚類在游動時胸鰭的擺動總呈現出一定的間歇運動，即胸鰭除了向后拍動，借助水的作用力推動魚體游動的加速階段之外，還包括回擺的過程以及2個擺動周期之間胸鰭相對靜止的過程，由于胸鰭并未提供向前的推力，魚體本身靠慣性作用向前滑行，因此以上2個階段為機器魚的滑行階段。在“加速—滑行”行為的研究中，萬宏[22]分析了魚尾鰭波動時的“加速1滑行”游動模式，發現在2個運動周期之間加入空白的滑行階段可提高機器魚游動時的效率，降低能耗比。將該模式代入胸鰭的運動中，可防止胸鰭擺動時動力階段產生的渦過早地被恢復階段胸鰭的回擺破壞掉，充分利用胸鰭后擺所產生的推力。

占空比，即在胸鰭的完整擺動周期中擺動階段需要的時間所占用的比例，其表達式如下：

(2)

式中：Tb為一個周期內胸鰭擺動階段的時間;Tc為滑行階段的時間。

為方便直觀的顯示占空比，現取2個周期內前后拍翼角的變化來表示胸鰭的運動，將占空比為1∶1、0.8∶1、0.6∶1、0.4∶1 時胸鰭的運動表示出來，如圖3所示，在周期一定的情況下通過改變胸鰭擺動的時間來調整占空比，在相同的時間內，不同占空比擁有相同周期。

2.2 魚體的計算模型

由于尾鰭采用了柔性波動，故運動規律的表示與胸鰭不同，設魚體的頂點處為慣性坐標系的原點，魚頭向魚尾的方向為X軸，垂直魚體水平方向為Y軸，垂直魚體豎直方向為Z軸，魚體的運動模型定義如下：

圖3 不同占空比下前后拍翼角的變化

(3)

式中：y為魚體在X位置的橫向位移，魚體的中軸線在一個周期內的波動曲線如圖4所示，由于魚體采用柔性波動，在UDF中運用節點運動將每個時間步內魚體中每一個節點的位移表示出來，并輸出到魚體模型上。

2.3 流場中的其他參數

機器魚在游動過程中的推力系數與側向力系數可分別由下式表達：

(4)

式(4)中：ρ為流體密度;Fx、Fy分別為魚體受到的推力與側向力;S為魚體的投影面積。

魚體的推進效率為平均輸出功率與平均輸入功率的比值，平均輸出功率為魚體運動速度與游動過程中平均推力的乘積，平均輸入功率為

(5)

式(5)中:σ為魚體的微元受到的力；V為該微元的速度，該乘積在魚體的表面積分。由此可得魚體的平均推進效率為

(6)

圖4 尾鰭波動曲線

3 流場的控制方程與數值方法

該流場的控制方程為非定常、連續、不可壓縮的三維Navier-Stokes方程，其表達式如下：

(7)

其中數值計算的網格如圖5所示，由于模型需要自推進來向前游動，設機器魚體長為L，計算域的尺寸為8L×1.5L×0.5L，且X方向的網格均勻分布，以保證在魚體游動過程中計算保持相同精度；流場網格劃分形式為四面體網格，模型的壁面采用三角形網格。

圖5 計算域的網格劃分

計算區域的邊界條件為：入口邊界與出口邊界的速度和壓力梯度均為0。由于在近壁面的流場中需要精度較高的結果，因此，本文實驗中所采取的湍流模型為k-ωSST模型，且采用低Re數校準，采用SIMPLEC算法對連續方程中的壓力和速度進行耦合，其中壓力、動量、湍流動能均采用二階迎風格式，擴散項采用二階中心差分格式，瞬態方程采用二階隱式。監測魚體在3個方向上的水動力系數與X方向上的速度與位移。

本文采用有限體積法(finite volume method,FVM),并利用動網格技術，通過剛體運動與節點運動來實現胸鰭的擺動和尾鰭的波動，并使用彈簧近似網格光順和局部網格重構方法以防止網格變形太大導致網格負體積的現象出現。該技術通過一套網格便實現了魚體、胸鰭運動的同時整體向前推進，不失精度的同時計算量相對較小。

魚體在計算域中采用牛頓定律以向前游動，魚體在t時刻前進的速度可由下式表示：

(8)

式中:Ab和Af分別為魚體與胸鰭的表面積;Fx(t)為機器魚在t時刻受到的X方向的合力，包括魚體收到的X方向的靜壓與左右胸鰭受到的力;m為機器魚的整體質量。

4 數值計算及結果分析

4.1 網格收斂性的驗證

在數值模擬中，為驗證網格的收斂性，取3種網格尺寸，對每一種網格進行同一條件的數值計算，最終計算出每一種網格的水動力系數如圖6所示，其中網格1、2、3分別對應尺寸粗、中、細，通過圖6可知，3種網格下的計算結果差異較小，因此，網格設置與數值計算方法合理且收斂，以下仿真結果均通過第二套中尺寸網格計算得到。

圖6 網格收斂性驗證

4.2 占空比對推力系數與速度的影響

通過對不同占空比的仿真計算，達到穩定巡游狀態下其中一個周期的機器魚推力系數與升力系數、其中一側胸鰭推力系數與所有周期的機器魚推進速度如圖7所示。

由圖7(a)可知，單個胸鰭的最大推力系數幅值隨占空比的減小而增大，占空比的改變導致了胸鰭擺動的頻率的改變，頻率越大，占空比越小，胸鰭推進效果越好，但是功率消耗也越大；由圖7(b)可知，機器魚的升力系數幅值同樣隨占空比減小而增大，且作用時間逐漸減少，在胸鰭滑行階段，機器魚的升力基本為0，可見機器魚的升力主要是胸鰭的擺動產生的；由圖7(c)可知，隨著占空比的減小，機器魚的最大推力系數逐漸增大，且在所有占空比下，推力系數正方向的數值略小于負方向，因此，在整體上機器魚的推力略大于阻力，機器魚的推力作用時間隨占空比的減小而減小，到滑行階段時，推力系數波動較小；由圖7(d)可知，隨著占空比的減小，穩定巡游速度逐漸增大，但是速度波動也逐漸增大，機器魚達到穩定巡游狀態的時間不隨占空比的改變而改變。

圖8為不同占空比下機器魚整個游動過程的平均推力、平均升力與升阻比。由圖8(a)可知，隨著占空比的增大，平均升力系數逐漸增大，且增大的幅度逐漸減小，占空比為1∶1即沒有中間滑行階段時平均升力系數最大，占空比0.8∶1時機器魚的平均升阻比最大；由圖8(b)可知，平均推力系數隨著占空比的增大而增大，結合圖7，占空比越大，推力系數峰值越小，但是平均推力系數越大。

圖7 不同占空比下推力系數和速度隨時間的變化

圖8 不同占空比下的平均推力系數、 升力系數與升阻比

4.3 流場速度矢量的分析

在胸鰭的運動規律中，取φRC=φFC=φFLC=0，取φR0=π/4、φF0=π/6、φFL0=π/12，取dφF=dφFL=π/6，尾鰭的波動方程中，取c0=0.02、c1=-0.08、c2=0.2，波數k=π、ω與胸鰭保持一致，在靜水中進行自推進數值模擬，對胸鰭擺動占空比為1∶1時通過在靜水中的數值模擬得到的速度矢量圖進行分析。

取某個周期的6個典型時刻，每個時刻的速度矢量如圖9所示。在t=T/6時，魚體尾鰭快速向上彎曲，在尾部左側產生了較大射流，從而產生一個較大的作用在魚體上的側向力，通過上個周期的運動，魚體側后方存在一個較大漩渦，該漩渦作用在魚體上產生了一定推力。隨后，尾鰭擺動幅度逐漸減小，尾鰭處的射流逐漸指向左后方，該射流的反作用力推動魚體向前游動，魚體受到的推力逐漸增大，漩渦脫落，開始遠離魚體；胸鰭開始展開，其周圍出現側向射流，到t=T/3時，胸鰭繼續通過擺動產生較小漩渦。

t=T/2時，尾鰭向左彎曲產生射流，該射流為魚體提供較大的側向力，而在魚尾后方，新的漩渦逐漸開始產生，增強魚體推力；胸鰭繼續向外側展開，得益于搖翼運動，胸鰭兩側壓差小，魚體受到的阻力較小。

t=2T/3至5T/6，尾鰭伸直，并開始向右側彎曲，側向射流進一步增強，導致魚體受到較大側向力，后方漩渦成型并逐步遠離魚體，產生射流，推動魚體游動；通過尾鰭作用新的漩渦產生。此時胸鰭向后劃水，產生一定推力。

t=T時，尾鰭和身體逐漸回到初始位置，射流再次在產生在尾鰭左側，后方的漩渦完全成型；胸鰭也回到初始位置，開始下一周期的運動，其后方存在的射流加強了推力。

圖9 流場的速度矢量

以上為一個周期內魚體周圍流場的速度矢量分析，通過尾鰭和胸鰭的擺動周期性的產生射流與漩渦，魚體間歇性的受到推力與側向力，推力推動了魚體向前游動，側向力總體上相互抵消，魚體不會發生明顯橫向運動。

4.4 流場壓力云圖的分析

圖10為上述周期同時刻機器魚周圍流場的壓力分布。t=T/6時魚體兩側壓差較大，且魚體后方有前一周期從魚體周圍脫落的低壓渦；t=T/3時，魚體兩側壓差進一步增大，且低壓區域向后運動，為魚體游動提供了推力，胸鰭兩側產生了不利于向前游動的壓差；t=T/2時，低壓區域進一步脫落，在魚體右側同樣產生了低壓區域；t=2T/3時，魚體左側的低壓區域脫落到魚體后方，魚體兩側壓差再次增大，為魚體提供了較大推力；t=5T/6時，胸鰭的擺動使兩側壓差變大，對魚體的前進產生了積極影響，魚體右側低壓區域開始脫落；t=T時，魚體左側再次產生低壓區域。

4.5 魚體周圍流場的三維渦結構

取上述同一周期同時刻魚體以及周圍流場，生成其三維尾渦結構，如圖11所示，三維渦結構的識別采用Q判據，渦結構的顏色反映了渦量的大小。

胸/尾鰭協同自推進機器魚在游動過程產生了十分復雜的流場渦結構，下面就穩定巡游狀態下一個周期的6個時刻做簡要分析。胸鰭和尾鰭的拍動產生了大小、方向不同的垂直于XY平面的渦環，這些渦環導致了機器魚的前游。

圖10 流場的壓力分布

t=T/6時，在尾鰭處，之前的運動產生了兩列渦環，尾鰭附近的渦強度較高，該渦為該時刻魚體收到的推力與側向力的主要來源，上個周期脫落的渦環逐漸遠離魚體，對魚體的作用逐漸減弱；在胸鰭的運動中，上個周期產生的渦環仍在產生作用，胸鰭的擺動正在產生新的渦環。由于尾鰭的擺動，在t=T/3時渦環逐漸增大，渦強度也逐漸增大，且主要集中在尾鰭的一側，一定程度上增強了尾鰭的側向力；在胸鰭處，上個周期產生的渦環完全脫落，新的渦環開始產生。

通過連續作用，t=T/2時，尾鰭處逐漸產生新的渦環VB1，該渦環逐漸增大并向機器魚的左側緩慢移動；兩側胸鰭處也產生了新的垂直于XY平面的渦環VF1與VF2，此時胸鰭與尾鰭產生的渦分別作用在機器魚上，暫未發生渦的融合現象。隨后VB1向后方移動，到t=2T/3時，尾鰭處的新渦VB2誕生，VB1的頭部與VB2的尾部相連接，且VB1與VB2在空間上相互垂直，共同對魚體產生作用；胸鰭處VF1與VF2同時開始從胸鰭前緣與已經脫落的渦環一樣，慢慢向尾鰭方向移動，逐步靠近尾鰭處的渦。

至t=5T/6，VB2逐漸增大，胸鰭的第一個渦環逐漸與VB1融合，魚體受到的推力被加強。隨著VB2向側后方移動，t=T時，右側胸鰭的渦開始與VB2融合，VB1與VB2開始逐漸分離并脫離魚體，并向后移動，尾鰭通過運動開始產生新的渦環，胸鰭處，隨著胸鰭的擺動，VF1與VF2強度增大，新的周期開始。

圖11 流場的三維渦結構

4.6 實體實驗驗證

為驗證數值計算的準確性，通過實驗室現有機器魚進行了靜水中的直游實驗，并通過超聲波傳感器監測機器魚的前進方向實時位移數據，實驗過程共持續8 s，其中6個時刻機器魚位置及運動軌跡如圖12所示，機器魚在流場中前進的同時，由于尾鰭的周期性波動，產生較小的橫向位移，但總體上呈現直游行為。CFD仿真與實體實驗中機器魚前進位移隨時間的變化如圖13所示，相同時間機器魚的實際位移均小于仿真結果，其主要原因為實驗中機器魚外殼的近親水性材料、實驗時的壁面效應以及外界因素的干擾增加了游動阻力；可在機器魚外殼填充疏水性材料、增加流場尺寸以及優化實驗流程來減阻，以改進實驗方法。實驗總體誤差在期望范圍內，且總體游動趨勢一致。

圖12 機器魚游動過程

圖13 機器魚位移隨時間的變化

5 結論

采用CFD方法數值模擬了不同占空比下仿生機器魚在水平面內的自主推進運動，其中仿生機器魚帶有一對分離式的3自由度胸鰭，其中胸鰭采用剛性擺動，機器魚體采用柔性波動。通過分析不同占空比機器魚自主推進過程中的魚體游動性能與機理、水動力性能、流場的速度矢量、壓力云圖與三維渦結構，得出了以下結論：

1) 在機器魚的自主推進游動過程中，通過胸尾鰭的協同推進，魚體胸鰭后方周期性的產生射流與漩渦，通過射流與漩渦的共同作用，為機器魚提供了向前游動的推力。

2) 在機器魚的整個擺動過程中，胸鰭和尾鰭兩側均有較大壓差，高壓區域脫落到魚尾附近推動機器魚向前游動。

3) 在每個周期的游動中，流場的渦結構都產生了復雜的變化，魚尾的擺動產生了卡門渦街，并且胸鰭產生的渦向后脫落，逐漸與尾鰭的渦融合，增加了渦的強度。

4) 對比不同占空比下機器魚整體和胸鰭產生的推力系數，發現隨著占空比的減小，推力系數的振幅逐漸增大，但是持續時間相應的減少，到滑行階段時，推力系數的波動較小。在速度的對比中，發現隨著占空比的減小，機器魚的穩定巡游速度增大，但是速度波動較大，占空比越小，魚體的推進效果越好，但是占空比的減小也導致了功率消耗的增加。