二維多關節機器魚遠場流速和攻角預測研究*

劉科顯，孟慶浩，徐雪寒

(天津大學電氣自動化與信息工程學院，天津市過程檢測與控制重點實驗室，機器人與自主系統研究所，天津 300072)

流場環境感知是水下機器人基本功能，而人造側線系統是目前這一領域研究的熱點之一。側線是一種魚類獨有的器官，可以感知流場壓力、流速、溫度和鹽度等流場信息[1]，也可以用來辨別附近是否存在障礙物。

受側線的結構和功能啟發，研究人員已提出多種人造側線系統。Sharif等[2]利用離子聚合物與金屬復合材料研制一種可以覆蓋在機體表面的微型壓力傳感器；Abels等[3]設計了一種基于壓力應變片的機體表面微型流速傳感器，可較為精確地感知流速；李淵等[4]設計了能有效感知水流信息的“人工側線”系統，李振等[5－6]設計一種以魚類側線為基礎的MEMS水聽器，并采用有限元技術對封裝結構進行優化設計。

在應用方面，吳乃龍等[7]詳細分析了魚類側線功能，并實現了水流場識別；Zheng等[8]利用壓力傳感器捕捉流場中的渦街信息，實現了對相鄰仿魚機器人的位置和運動狀態的感知；Xu等[9]采用最優權重分析算法優化了“人工側線”系統的傳感器排布方式；Carryon等[10]將感知反饋策略應用于波動推進驅動器，提高了仿魚機器人的運動性能。劉貴杰團隊[11－13]建立了基于壓力傳感器的靜態“人工側線”系統，采用機器學習技術實現流場流速、物體位置等信息的流場結構反演。林興華[14]等采用CFD方法獲得不同類型障礙物的流場特征，采用支持向量機預測不同類型的流場目標。Salum?e等[15]采用人造側線系統辨識機體運動速度等運動狀態信息。Tallapragada等[16]分析了水翼擺動推進工況下翼型周圍流場分布狀態，可進一步用于分析仿生機器魚周圍流場規律。目前，對于魚類推進機理的研究主要以鰻鱺和鲹科魚類為主，但一般在單剛體或者柔性結構上安裝人造側線系統。

為擴展人造側線系統的應用范圍，本文以仿鲹科魚類游動的二維多關節機器魚為研究對象，將壓力傳感器陣列人造側線系統應用到具有多關節特征的機器魚結構上，研究在廣闊水域中多關節機器魚對遠場流速和攻角的感知方法。以Joukowski翼型作為二維機器魚包絡外形，在機器魚的頭部和魚體兩側設置11個壓力監測點，通過組合各壓力監測點輸出值給出了四種測量模型。應用勢流理論和伯努利方程分析Joukowski翼型表面壓力和流速，并結合Fluent軟件分析遠場流速和攻角與四種測量模型輸出值之間的關系，采用最小二乘多項式回歸方法建立四種測量模型的回歸方程。最后采用SIR粒子濾波實現以四種測量模型輸出值對遠場流速和攻角的預測。

1 擺動狀態Joukowski對稱翼型體表壓力分析

設翼型處于無限大的流場中，遠場流速為V∞，翼型攻角為α。定義流場參數為Λ＝(V∞，α)。通過式(1)所示的保角映射，轉化到復平面內的圓盤分析，即

式中，z表示翼型上某點，ξ為z在復平面上的映射，b0表示復平面圓盤圓心位置，此處可定義成b0＝R－｜ξ0｜，R表示轉換后圓盤直徑，ξ0表示圓盤圓心的位置，故Joukowski對稱翼型在復平面內可以用R和ξ0表示。則翼型周圍流場復勢可以表示為[16]:

式中，W1(ξ)為翼型固定時周圍流場的復勢，定義Γ為速度環量，根據勢流理論和Kutta條件可以寫出:

式(2)中:設Ω為翼型擺動的頻率角速度，WΩ(ξ)表示與翼型擺動相關的復勢，該項可表示為[16]:

式(2)中:W(k)v(ξ)為第k個渦的速度勢，根據Milne-Thomson圓定理，該項可表示為[16]:

式中，Γk為第k個渦的速度環量，ξk表示第k個渦的位置。翼型表面共軛速度場是速度勢的位置梯度，設翼型上某點流速為v(z)＝u＋iv，則共軛流速ˉv2(z)在復平面內可表示為:

式中，ξ(z)為復平面內一點到z平面的逆映射，可表示為

假設翼型擺動的角速度為Ω＝2πf，代入式(6)中并簡化可得

式中，H1(z)、H2(z)和H3(z)為關于翼型上某點位置z的函數，H4(z)是表征渦位置的參數。將上式中的指數項泰勒展開并取前兩階，則式(8)所述的ˉv2(z)可表示為

式中，ˉv3(z)的模為位置z處的流速，將式(9)代入伯努利方程，則翼型上任意一點的靜壓可表示為

式中，p(z；Λ)為流場參數Λ下翼型上位置z處的體表壓力；C為流場常量；ρ為流體密度。

為了消除流場常量C的影響，設翼型上任意指定兩點z1和z2的壓力差為:

從式(10)中可以得出，翼型上任一點壓力可以分成兩部分，一部分是由流場參數Λ決定，第二部分是由流場參數Λ和擺動頻率f耦合作用。

多關節機器魚周圍流場較為復雜，難以根據勢流理論建立體表壓力模型，需采用商業CFD軟件Fluent分析多關節機器魚體表壓力與遠場流速V∞和攻角α特性。

2 二維多關節機器魚仿真模型及魚體波動方程

2.1 二維多關節機器魚幾何模型

本文所采用的二維多關節機器魚模型如圖1(a)所示，其中B0為魚頭，B1、B2和B3為魚身，B4為機器魚尾鰭。魚頭和魚身為Joukowski對稱翼型，尾鰭為NACA0013翼型，魚體總長為320 mm，其中魚身長度168 mm，在魚體周邊均勻分布11個壓力監測點。迎著魚頭觀察，左側壓力監測點分別命名為s11、s21、…、s51，右側壓力監測點分別命名為s12、s22、…、s52，在魚頭正中處設置壓力監測點s0。

采用如圖1(b)所示組合網格方案，機器魚周圍流場采用三角形網格，遠處流場采用四邊形網格。

圖1 多關節機器魚模型

2.2 魚身波動方程

魚身部分波動規律需符合Lighthill[17]提出的鲹科魚類魚體波方程擬合，如式(12)所示。

式中，ybody為魚體的橫向位移，x為魚體的縱向位移，c1和c2為波幅包絡線系數，k＝2π/λ為魚體波波數，ω＝2πf為魚體波角頻率。其中λ表示魚體波波長，f為魚體波擺動頻率。

將魚體波方程在一個周期內離散，可得到如式(13)所示的一組離散后的魚體波曲線簇。

式中，M為魚體波在一個周期內的離散數量，本文取M＝18，i為曲線簇中的曲線編號，取值范圍為0～M－1，定義c1＝0.12，c2＝0.21，波長λ取1.5倍體長，假設構成魚身的各段旋轉中心位于魚體波曲線上，則可得式(14)表示的各連桿轉動規律，

式中，φ1max、φ2max、φ3max和φ4max分別為魚身B1、B2和B3段以及尾鰭B4最大轉角；φ1為B2相對于B1滯后的相位差；φ2為B3相對于B2的相位差；f為擺動頻率；t為擺動時間；p為指數系數；Lbody為魚身長度。

3 多關節機器魚體表壓力分析

3.1 壓力系數測量模型

各壓力監測點數值隨著魚體波動呈周期性變化，需對采樣數據做式(15)所示的時間均值處理以及無量綱化處理。

式中，p(si；Λ)為經過時間均值處理后的壓力值，Cp(si；Λ)為監測點si在流場參數為Λ條件下測量的無量綱壓力系數；p∞為遠場壓力或參考壓力，此處選一個大氣壓；U*為參考速度，數值上等于魚體波波長λ。建立以下四種測量模型，分別命名(1)單點傳感器，標記s；(2)交叉傳感器，標記為sopp；(3)相鄰傳感器，標記為sadj，(4)對稱傳感器，標記為ssym。按照式(15)所給出的無量綱壓力系數Cp(si；Λ)的定義，四種測量模型各測點輸出值可以寫成式(16)的形式。

式中，s、sopp、sadj和ssym分別表示單點傳感器、交叉傳感器、相鄰傳感器和對稱傳感器模型；βi(i＝1，…，32)為四種測量模型測點輸出值。

3.2 靜水環境下機器魚體表壓力特征

設二維機器魚處于靜水環境，魚體按照式(14)擺動，設擺動頻率為1 Hz，機器魚一個擺動周期內壓力場云圖以及體表壓力時空分布情況分別如圖2和圖3所示。

圖2 單個擺動周期內的壓力場云圖

圖3 各測量模型時空分布

圖2中隨著云圖顏色越深，流場壓力越大，可看出，尾鰭附近存在高壓區，魚體同一側從尾部到頭部壓力逐漸減小，魚體兩側流場交替出現高壓區和低壓區，具有時空對稱性。

從圖3(a)和3(b)中可以看出單點傳感器得到的壓力系數具備時空對稱性。從圖3(c)和3(d)中可以看出，去掉流場常量的影響后，相鄰傳感器和對稱傳感器相比于單點傳感器對壓力系數波動有一定的抑制作用，圖3(c)所示的對稱傳感器所得到的壓力系數在空間和時間上的分布更為均勻。

圖4為兩個擺動周期部分單點傳感器和相鄰傳感器測量值隨時間變化曲線。壓力系數周期與魚體波動周期相同。圖4(b)中關于魚體中軸線對稱的一對單點傳感器模型測得的壓力系數相位差相當于半個擺動周期，魚體同側的單點傳感器相位差基本為零。圖4(c)中單點傳感器壓力系數波動幅值從魚體頭部到魚體尾部逐漸增大。圖4(d)所示的相鄰傳感器壓力系數時域曲線波動幅值比單點傳感器低一倍左右。

圖4 各測量模型測量值隨時間變化曲線

3.3 不同流場流速機器魚體表壓力特征

設遠場流速范圍為0到2 m/s，魚體擺動頻率為1 Hz，魚首攻角為2°，計算不同流速流場條件下各測量模型監測的壓力系數，結果如圖5所示。

各測量模型測得的壓力系數Cp基本是遠場流速V∞的二次函數，在指定的流場參數范圍內具備單調性，圖5(a)、5(b)和5(c)所示的單點傳感器測得的壓力系數在流速為0.6 m/s到1.2 m/s區間內有明顯波動。本文分析采用“網格重構”法構建動網格技術，重構時造成數值計算偽擴散，從而形成魚體表面壓力的波動，可以視作誤差處理。

圖5 壓力系數Cp與來流速度變化曲線

由于消除一定誤差以及流場常數的影響，圖5(d)、5(e)和5(f)所示的組合傳感器模型測得的壓力系數相比于單點傳感器更為平穩。

為了評價四種測量模型各測點測量值離散程度，繪制圖6所示的四種測量模型各測點壓力系數標準差隨遠場流速V∞變化趨勢。隨著遠場流速V∞增大，單點傳感器和交叉傳感器各測點測量值標準差迅速增大，相鄰傳感器和對稱傳感器測量值相對平緩，更適合處理遠場流速變化劇烈的流場環境。

圖6 不同流速下壓力系數的標準差

3.4 不同攻角機器魚體表壓力特征

設遠場流速V∞＝1 m/s，攻角范圍為－16°到16°，魚體擺動頻率為1 Hz，各測量模型得到的壓力系數如圖7所示。

從圖7的曲線趨勢上看，壓力系數Cp和攻角α之間呈高于兩次的函數關系。圖7(a)和圖7(b)所示的相鄰傳感器和單點傳感器在指定攻角范圍內不具備單調性，其中后者在攻角為0兩側單調性相反，圖7(c)和圖7(d)交叉傳感器和對稱傳感器測量模型在全域上具備單調性。

圖7 壓力系數Cp與攻角關系曲線

從數值波動的角度來看，圖7(a)、7(c)和7(d)所示的相鄰傳感器、交叉傳感器以及對稱傳感器測量數據比單點傳感器平滑，交叉傳感器和對稱傳感器可以利用魚體兩側信息從而判斷攻角的正負。

根據圖8所示的不同測量模型測量值標準差隨攻角變化曲線來看，無論是迎流側還是背流測，相鄰傳感器和單點傳感器標準差較小，對稱和交叉傳感器測量值標準差隨著攻角增加迅速增大。對于攻角波動較大的情況，例如機器魚做復雜曲線運動，此時可以用交叉或對稱傳感器預測攻角正負，用相鄰傳感器測量具體數值。

圖8 不同攻角下壓力系數的標準差

3.5 測量模型輸出值回歸方程

根據上述分析，魚體擺動頻率一定的情況下，魚體體表壓力是遠場流速V∞和攻角α的函數。在CFD計算的基礎上擬合來流場參數及波動參數與體表壓力系數之間關系，設壓力系數函數為Cp(V∞，α；i)。

根據第3.1和3.2節分析，以式(10)、式(11)和式(15)為基礎可知，在擺動頻率f已知時，體表壓力系數可近似視為關于攻角α四次和遠場流速V∞二次的函數。故本節根據上述計算得到的數據，基于最小二乘法，以V∞和α為變量構建壓力系數Cp(V∞，α；i)多項式回歸方程。

由于觀測點關于魚體中軸線對稱分布，對于單點傳感器和相鄰傳感器只需構建迎流側壓力系數函數即可。設基準擺動頻率為1 Hz，以對稱傳感器ssym3、交叉傳感器sopp32以及迎流側的單點傳感器s32和相鄰傳感器sadj22為例，構建Cp(V∞，α；i)回歸模型，結果如圖9所示。

圖9 四種測量模型擬合結果

從表1中可以看出采用以式(10)和式(11)形式構建的多項式回歸方程可以較為精確地擬合各采樣點，相鄰傳感器擬合效果要優于其他測量模型。

表1 各測量模型擬合度

4 基于SIR粒子濾波的流場參數預測

將3.5節提出的回歸方程作為觀測方程，估計流場參數，即x＝(V∞，α)。設實際測量值為yi＝βi＋ηi，其中ηi～N(0，σ2i)是第i個測量模型的噪聲，呈均值為0，方差為σ2i的高斯分布，不同測量模型之間的高斯噪聲互相獨立，定義y＝(y1，y2，…yN)為各測量模型所測量值的集合，則流場參數x的后驗概率為

式中，κ為歸一化因子，p(y｜x)為似然函數，p(x｜y0)為先驗概率分布，設為均勻分布。

為了融合多個測量模型所獲得的信息，采用多元高斯分布，對于單個測量模型其似然函數為

由于各測量模型高斯噪聲互相獨立，因此流場參數x的后驗概率密度為

則系統的預測方程為

系統的更新方程為:

式中，Y(t)＝{y(t)，y(t－Δt)，…，y(0)}，設先驗概率的初始條件p(x(0)｜Y(－Δt))為高斯分布。

設遠場流速和攻角呈正弦變化，其中遠場流速波動幅值為0.5 m/s，攻角波動幅值為1.5°，周期均為2 s，即:

觀測方程為:

根據式(20)和式(21)定義的預測模型，應用SIR粒子濾波預測流場參數，即遠場流速V∞和攻角α，測量模型各測點獲取的壓力系數為觀測量，通過數值計算得出，遠場流速V∞由入口速度給定，攻角α由機器魚自身轉動設定，粒子濾波相關參數見表2。

表2 粒子濾波參數

為去除測點位置帶來的影響，選取s32、sopp32、ssym3和sadj32四個同一位置處的各測量模型輸出的回歸方程作為觀測方程，重復試驗50組，設機器魚擺動頻率為1 Hz，仿真時長6 s，仿真結果和對應的偏差曲線如圖10和圖11所示。

圖10 攻角粒子濾波預測效果和偏差曲線

圖11 遠場流速粒子濾波預測效果和偏差曲線

從圖10和圖11中可以看出，四種測量模型輸出結果均可以較為準確地預測流場參數，預測結果基本跟隨真實狀態。根據偏差圖，在攻角α的估計上四種測量模型的偏差基本一致，對于V∞的估計上相鄰傳感器偏差最小，交叉傳感器和對稱傳感器次之，單點傳感器最大。由于相鄰傳感器、交叉傳感器和對稱傳感器三種測量模型考慮魚體表面兩個位置的壓力系數，一定程度上抵消了部分流場噪聲。

為評估多傳感器組合應用對流場參數的估計效果，以相鄰傳感器sadj12、sadj22，sadj32和sadj42為例組合預測。為反映位置因素對觀測產生的影響，從魚首到魚尾依次加大測量噪聲，仿真計算50次，預測效果和偏差曲線如圖12所示。計算多傳感器組合應用與相鄰sadj32傳感器預測偏差的絕對平均偏差，結果如圖13所示。

將圖12與圖10(d)和圖11(d)對比，兩者在對遠場流速V∞的估計上精確程度相近，但較大提升了對來流角度預測精度。從圖13中可以看出兩個參數預測上多傳感器組合應用的穩定性較單一傳感器略有提升。

圖12 多傳感器組合預測流場參數

圖13 流場參數預測偏差的平均絕對偏差

5 結論

本文結論如下:①以勢流理論和伯努利方程為基礎分析得到Joukowski對稱翼型表面壓力與遠場流速、攻角和擺動頻率相關，并給出具體表達式；②在外包絡線為Joukowski對稱翼型的二維多關節機器魚基礎上建立四種測量模型并用CFD技術分析，各測量模型優缺點，建立回歸方程。綜合來看，相鄰傳感器模型測得的壓力系數效果較好；③以SIR粒子濾波為基礎的流場參數預測仿真表明所建立的測量模型可以實現對遠場流速和攻角的預測。多個測量模型融合后預測效果的平均絕對偏差小于單個測量模型預測結果。

本文提出的方法可以預測較為寬闊的區域的遠場流速和攻角，為多關節機器魚的流場適應性控制提供信息，也為多關節機器魚探測復雜流場環境方法提供解決思路。