張開升 鄒慶彪 趙波 張保成

(中國海洋大學工程學院，山東青島266100)

仿生機器魚作為水下機器人的一種，具有高效、高速和靈活的特點，逐漸成為探索水下世界的重要手段。與傳統(tǒng)的多剛體仿生機器魚相比，柔性仿生機器魚由柔性材料制造而成，具有更好的柔順性、魯棒性、高效性和低擾動性，能夠更好地模擬魚類的游動姿態(tài)和適應更加復雜多變的流場環(huán)境，逐漸成為國內外研究的熱點。但當前國內外對柔性仿生機器魚的研究還主要集中于采用實驗手段進行柔性仿生機器魚的功能實現(xiàn)研究，如Alvarado［1］制作了基于硅膠材料的柔性機器魚，Marchese等［2］制作了由流體系統(tǒng)驅動的柔性機器魚，Wang等［3］設計制作了四關節(jié)機器魚，并致力于實現(xiàn)其最大游泳速度和更高的推進效率，而在魚體設計參數(shù)對魚體游動性能的影響規(guī)律方面研究還存在不足，因而也導致仿生機器魚樣機與真實魚類的游動性能存在較大差異［4］。另外，通過實驗的方法對柔性機器魚進行研究存在研究周期長、研究成本高的缺點，且在樣機制作完成后其設計參數(shù)不能改變，從而導致無法開展較為全面的研究，致使研究效率較低。因此，通過數(shù)值仿真的方法對柔性仿生機器魚進行研究是一種可行且高效的研究手段。目前，許多學者對仿生機器魚進行數(shù)值仿真研究，如Xia等［5］基于CFD(計算流體動力學)，采用UDF(用戶自定義函數(shù))和動網(wǎng)格技術對仿生機器魚自主推進直線運動進行研究;Kern等［6］建立了鰻魚的三維自主游動模型，并研究了魚體周圍受力的分布情況和尾跡渦形成與脫落的具體過程;楊亮［7］采用CFD的方法模擬仿金槍魚尾鰭擺動時產(chǎn)生的尾渦特點，并分析其與游動性能之間的關系;胡亞南［8］以浸入邊界法為理論基礎，通過自主開發(fā)的算法程序研究仿生魚的游動機理;劉鵬等［9］采用CFD的方法研究了柔性翼推力及效率;李寧宇等［10］基于浸入邊界有限體積法研究拍動翼運動學、渦動力學和力的產(chǎn)生之間的關系;張開升等［11］采用CFD的方法研究了仿生機器魚的多機體協(xié)同推進效率。但上述研究并沒有考慮因為流固耦合作用而引起的魚體變形，因此造成仿真結果存在較大誤差。

1964年，Noh［12］提出的耦合歐拉 -拉格朗日(簡稱CEL)法能分析物體運動導致的流體大變形問題，在解決流體與結構之間的復雜相互作用問題時具有強大優(yōu)勢，如王永虎等［13］基于CEL法研究某型飛機撞入水面的流固耦合過程;王懿等［14］基于CEL法研究船舶拋錨入泥的流固耦合過程。因此，CEL法能夠分析柔性仿生機器魚在游動過程中魚體的柔性變形以及魚體與流體之間的流固耦合作用，能夠為柔性仿生機器魚的研究提供一種新的研究思路。本文采用CEL法建立柔性仿生機器魚自主游動的流固耦合模型，研究魚體剛度、擺動幅值、擺動頻率對柔性仿生機器魚巡游速度的影響規(guī)律，據(jù)此得到仿生機器魚達到較優(yōu)巡游速度的設計參數(shù)組合。

1 耦合歐拉-拉格朗日法

傳統(tǒng)的分析方法，如拉格朗日法 (Lagrangian)，可以用來分析物體的受力及變形，但網(wǎng)格須與分析目標的變形保持一致，即材料與單元的變形必須保持一致。該方法可以簡化求解，在小變形的情況下能夠得到較為準確的節(jié)點及位移等信息，但變形較大時會導致網(wǎng)格畸形，從而無法收斂致被迫中止計算。而歐拉分析法 (Eulerian)中網(wǎng)格的位置屬性、形狀屬性不會變化，包括流體在內的材料可在網(wǎng)格中流動，不會出現(xiàn)無法收斂的問題，因此該方法適于分析和處理變形較大的問題。但其網(wǎng)格的固定致使無法準確地捕捉物體的邊界信息，且不能得到結構的準確受力情況。兩種算法中的網(wǎng)格單元特征及連續(xù)介質變形情況如圖1所示。

圖1 網(wǎng)格特征及連續(xù)介質變形情況Fig.1 Mesh characteristics and continuum deformation

1.1 相互作用及控制方程

歐拉材料的質量守恒方程、動量守恒方程、能量守恒方程及連續(xù)性方程在歐拉坐標系下進行時間積分，且材料之間的相互作用使用黏性行為，拉伸應力可以在兩個歐拉材料之間穿過界面進行傳遞，并且在這些界面上不會發(fā)生滑動。

CEL法的接觸公式基于增強浸入邊界法，為了模擬物體之間的相互作用，兩套網(wǎng)格之間的接觸在有限元軟件Abaqus中使用基于罰函數(shù)方法的通用接觸。仿真分析時，系統(tǒng)能自動地計算并且追蹤拉格朗日結構和歐拉材料之間的界面，而不需要為歐拉域生成一個協(xié)調網(wǎng)格。該方法中接觸面間的法向接觸壓力與兩個接觸面網(wǎng)格節(jié)點之間的法向間距h有關，如圖2所示。

圖2 CEL法的接觸判定示意圖Fig.2 Contact determination schematic diagram of CEL method

如果發(fā)生節(jié)點對物質界面的貫穿，界面力就會傳遞到歐拉材料的錨定點上，接觸壓力如下所示:

式中:Fp為接觸壓力;K為罰剛度系數(shù)，取決于拉格朗日和歐拉材料的特性;h為錨點與網(wǎng)格種子之間的距離。

1.2 顯式動力分析

在有限元軟件Abaqus中使用CEL法進行仿真計算，采用的是Abaqus/Explicit顯式分析求解器。該求解器采用顯式中心差分時間積分法則［12］，可用于分析柔性魚在流場中的擺動行為。該分析方法的穩(wěn)定性限制了Abaqus/Explicit求解器所能采用的最大時間步長，因此為了保證求解的精度和可靠性，該求解器默認使用自動時間增量，使時間增量盡可能地接近而且又不超過穩(wěn)定性限制。時間增量與整個模型的固有頻率相關。

1.3 材料屬性的定義方式

流體的材料屬性通過狀態(tài)方程 (簡稱EOS)中的Us-Up狀態(tài)方程，來描述模擬流體動力學材料或幾乎不可壓縮流體和描述材料的沖擊力、體積和能量特性［15］。在CEL法中，流體的材料被描述為黏性可壓縮牛頓流體，在黏性的作用下流場會產(chǎn)生渦流，可用下式描述:

式中:σ為Cauchy應力張量;p為壓強;η為剪切黏度;˙e為應變率。

2 基于CEL法的流固耦合有限元模型

在Alvarado［1］提出的鮪科和鲹科魚類的外形方程基礎上，利用三維軟件建立仿生機器魚的三維模型，其長度為0.5 m;整個魚體為拉格朗日體。由于本文只考慮柔性仿生機器魚的直線巡游運動，魚頭擺動的幅度及變形較小，且通常頭部的剛度較大，因此，為簡化模型，將仿生機器魚有限元模型頭部設置為剛體并約束其只能沿X方向運動;由于柔性尾鰭的運動模式通常是隨著魚體的運動而被動擺動，因此將尾鰭設置為剛體并使尾鰭繞其與后頸部的連接處做擺動和平動，如圖3所示。

圖3 仿生機器魚三維模型示意圖Fig.3 Schematic diagram of three-dimensional model of bionicfish

魚體材料始終在流場壓力和魚體運動的作用下變形，尾鰭始終隨著尾部被動擺動，且被動擺動的幅度與魚體材料的變形和流場有關，最終仿生魚在流固耦合的作用下完成最終的游動表現(xiàn)。

流體域必須選擇為歐拉體。根據(jù)仿生機器魚的尺寸和可能的游動范圍，將流體域尺寸設置為:長2.50m、寬0.35m、高0.25m。

柔性魚體的材料屬性主要考慮密度ρ、彈性模量E。根據(jù)文獻 ［16］，流體域中水被定義為近似不可壓縮的黏性流體，材料的基本屬性參數(shù)如下:線性沖擊波波速Us的系數(shù)c0=1485m/s，粒子速度Up的系數(shù)s=0，材料常數(shù)Γ0=0，密度ρ=1000kg/m3，黏度 μ =1 ×10-3Pa·s。

為盡可能地避免魚體材料由于彈性模量較小而在流體中發(fā)生拉伸變形，在仿生機器魚沿魚體方向的中間截面上添加一層具有較大彈性模量的膜單元層，如圖4所示。

圖4 仿生機器魚模型中的薄膜層Fig.4 Membrane layer in biomimetic fish model

使用波動方程來描述魚類的波狀運動，因此在有限元軟件Abaqus中，選擇仿生機器魚體干上的若干點，根據(jù)波動方程獲得每個點的側向位移數(shù)據(jù)并輸入軟件中。波動方程建立的基礎是仿生魚做直線連續(xù)運動，尾部做被動擺動的運動，該運動模型有如下假設:①橫向行進波在魚體每一個截面上都是一維的，沿著頭部向尾部傳播，且沒有沿行進方向的扭動;②在波動游動中，通過身體伸縮來實現(xiàn)這種運動［17］。根據(jù)波動方程可知，魚身中部的擺動幅度遠小于尾部，因此在尾部選取更多的點，每個點之間的運動幅度通過樣條差值獲得。

根據(jù)筆者前期的研究［11］，基于牛頓歐拉法和遺傳算法獲得了本文有限元模型中使用的基本波動方程:

式中，h(x，t)為從魚頭到魚尾方向體干上坐標為x處t時刻的側向位移。

由CEL法可知，拉格朗日體的表面自動和歐拉體進行接觸，通過對拉格朗日體表面進行積分得到節(jié)點力 (歐拉體施加的壓力)。拉格朗日單元相當于給歐拉單元施加了流動約束，部分拉格朗日體通過覆蓋歐拉網(wǎng)格的體積和面，將流動約束施加給歐拉物質，因此歐拉單元的尺寸和拉格朗日單元的尺寸之間的比例對耦合效果具有一定影響，即歐拉單元的網(wǎng)格質量對計算結果存在影響。除此之外，罰剛度系數(shù)也會對計算結果產(chǎn)生影響。根據(jù)現(xiàn)有的硬件條件，網(wǎng)格收斂性分析對比組的設置及結果如表1所示。由表1可知，網(wǎng)格組5、組4的計算結果趨向收斂狀態(tài)。為減少計算成本以提高效率，本文采用組2的方案，因此也導致最終計算的結果與收斂值相比略有偏差。

表1 網(wǎng)格收斂性分析結果對比Table 1 Comparisons of mesh convergence analysis results

有限元模型的網(wǎng)格數(shù)量及網(wǎng)格類型見表2。

表2 魚體及流體域網(wǎng)格數(shù)量及類型Table 2 Number and type of mesh of fish body and fluid domain

柔性仿生機器魚自主游動的流固耦合有限元模型如圖5所示。

圖5 柔性仿生機器魚自主游動的流固耦合有限元模型Fig.5 Fluid-Structure-Interaction finite element model for selfpropulsion flexible bionic fish

3 設計參數(shù)對仿生機器魚巡游速度的影響

3.1 設計參數(shù)組合的確定

本文采用正交試驗來確定設計參數(shù)組合。正交試驗能夠以較少的試驗成本得到最優(yōu)的試驗方案和結果，在式 (3)的基礎上，本文的正交試驗主要考慮魚體剛度、擺動幅值和擺動頻率這3個參數(shù)。通過密度、彈性模量和泊松比定義魚體剛度，彈性模量定為10～50MPa。由鲹科魚類的擺動幅值而將擺動幅值定為0.10l(l為體長，0.5 m)至0.20l。擺動頻率定為1.72～5.16 Hz。將每個參數(shù)分為5個水平 (見表3)。

表3 因素水平表Table 3 Factor-level table

對三參數(shù)五水平選用正交表L25(56)確定不同參數(shù)組合和試驗安排，共25組試驗。

3.2 各個設計參數(shù)對巡游速度的影響

由正交表，利用有限元軟件Abaqus進行仿真計算，使用極差法對仿真結果進行統(tǒng)計，并得到各設計參數(shù)對巡游速度的影響規(guī)律。

3.2.1 魚體剛度對巡游速度的影響

將正交實驗的結果擬合可得到各個設計參數(shù)對游動性能的影響規(guī)律。魚體剛度對巡游速度的影響規(guī)律如圖6所示。由圖6可知，巡游速度呈現(xiàn)先增后減，然后趨于平穩(wěn)的趨勢，最大值出現(xiàn)在30MPa附近，在40至50MPa之間穩(wěn)態(tài)巡游速度呈現(xiàn)略微上升趨勢。

圖6 仿真樣本內魚體剛度對巡游速度的影響Fig.6 Effect of fish stiffness on speed of cruise in simulated samples

3.2.2 魚體擺動幅值對巡游速度的影響

魚體擺動幅值對巡游速度的影響規(guī)律如圖7所示。由圖可知，巡游速度總體呈現(xiàn)先增后減的趨勢，在0.13l附近具有較高的巡游速度。

圖7 仿真樣本內擺動幅值對巡游速度的影響Fig.7 Effect of amplitude on speed of cruise in simulated samples

3.2.3 魚體擺動頻率對巡游速度的影響

魚體擺動頻率對巡游速度的影響規(guī)律如圖8所示。由圖8可知，巡游速度總體呈現(xiàn)隨擺動頻率單調遞增的趨勢，在較高的擺動頻率下可以獲得較大的推進速度。

圖8 仿真樣本內擺動頻率對巡游速度的影響Fig.8 Effect of swing frequency on speed of cruise in simulated samples

3.2.4 較優(yōu)設計參數(shù)組合

由圖6-8可知，當魚體剛度為30 MPa、魚體擺動幅值為0.13l、魚體擺動頻率為5.16 Hz時，柔性仿生機器魚能獲得較優(yōu)的巡游速度，且通過仿真結果可知該設計參數(shù)組合下仿生機器魚的巡游速度能夠達到約1.4l/s。

3.3 柔性魚體的變形

柔性魚體發(fā)生彎曲導致尾鰭的運動與魚體存在一定的滯后，從而影響兩者之間的協(xié)同擊水。對仿真進行觀測，發(fā)現(xiàn)在魚體剛度較小時，魚體后頸部出現(xiàn)壓縮彎曲現(xiàn)象，導致尾鰭運動的明顯滯后，且在擺動換向時尾鰭未完成當前方向的擺動即被魚體拖動向另一個方向運動，隨機選取一段游動過程如圖9所示 (t0為游動過程的某個時刻，T為擺動周期)。尾鰭與魚體不能良好地協(xié)同擊水，從而導致了性能指標較差。

圖9 魚體剛度較小時尾部出現(xiàn)的彎曲現(xiàn)象Fig.9 Bending phenomena in the tail with small stiffness

后頸部的彎曲程度會隨著魚體剛度的增加而逐漸減小，尾鰭運動滯后的程度也隨之逐漸減小。魚體剛度在30 MPa附近時，尾鰭與魚體能較好地協(xié)同運動，因而仿生機器魚能獲得較高的巡游速度。

為更清楚說明不同剛度對魚體彎曲的影響，給出不同魚體剛度在擺動幅值及擺動頻率分別為0.20l及5.16Hz時，相同時刻后頸部的彎曲程度，如圖10所示，其中左圖均為尾鰭擺動方向即將改變前的情況，右圖均為擺動方向改變后的情況。

圖10 尾部的彎曲情況Fig.10 Bending situation of fish tail

魚體擺動幅值和擺動頻率也會影響魚體的彎曲情況。對樣本內尾鰭擺動角度的最大值進行統(tǒng)計(見圖11)可以發(fā)現(xiàn):尾鰭最大擺動角度隨著擺動幅值的增加呈現(xiàn)上升趨勢，擺動幅值對尾鰭的最大擺動角度影響較為明顯，在擺動幅值為0.20l時，尾鰭最大擺動角度達到最大值，大于1.570 8 rad(90°)，造成尾鰭的位置在前進方向上超過尾柄，因此形成了明顯的負擊水行為，從而導致游動性能的下降;擺動頻率對尾鰭最大擺動角度的影響較小，隨著擺動頻率的增加，尾鰭最大擺動角度呈現(xiàn)下降的趨勢，在擺動頻率為1.72 Hz附近時達到最大值 (約1.38rad)，而在5.16Hz附近時達到最小值 (約1.32rad)。

圖11 仿真樣本內擺動幅值和擺動頻率對尾鰭最大擺動角度的影響Fig.11 Effects of amplitude and frequency on maximum swing angle of caudal fin in simulated sample

結合圖7、圖8、圖11可以看出，在最大擺動角度為1.32rad附近時，柔性仿生機器魚能達到較高的游動速度。

3.4 流場變化情況

選取新月形尾鰭中部的一部分，對經(jīng)過該部分邊緣的流線進行跟蹤，如圖12所示。由圖可知:在擺動過程中，尾鰭的擺動會使附近流體加速，從而逐漸形成閉合的環(huán)線，即逐漸形成渦，其旋轉速度和尾鰭擺動的角度和速度有關。

圖12 尾鰭中部形成的流線示意圖Fig.12 Schematic diagram of streamlines formed in the middle of caudal fin

仿生魚游動過程中，流體材料的某一瞬時壓力情況如圖13所示。由圖可知，尾鰭的擺動使流場產(chǎn)生高壓區(qū)和低壓區(qū)，尾柄處的流體壓力明顯較大，說明魚體通過擺動從流場獲得了較大的反作用力;尾鰭處的渦團逐漸脫落并散落在流場中，與另一側脫落的渦團相對應;魚體表面的流體壓力比尾鰭處的壓力小，前后的壓力差驅動仿生魚前進。

圖13 仿生魚游動過程的某一時刻Fig.13 A moment in the swimming process of bionic fish

4 結論

本文基于CEL法研究了柔性仿生魚的自主游動過程，首先建立了柔性仿生機器魚自主游動的流固耦合模型，然后在此基礎上，通過在有限元軟件Abaqus上進行仿真計算得到魚體設計參數(shù) (魚體剛度、魚體擺動頻率、魚體擺動幅值)對仿生機器魚巡游速度的影響規(guī)律，最后據(jù)此得到較優(yōu)的設計參數(shù)組合。在樣本范圍內，當柔性仿生機器魚具有一定的魚體剛度時能獲得較好的游動性能，對于本文的模型，當魚體剛度約為30 MPa時具有較好的游動性能;魚體剛度能影響魚體的彎曲程度，在此基礎上，魚體擺動幅值和擺動頻率能夠進一步影響魚體的彎曲程度和尾鰭的最大擺動角度，從而影響游動性能，對于本文的模型，魚體擺動幅值約為0.13l、擺動頻率約為5.16Hz時具有較好的游動性能。因此，對于本文的模型，基于CEL法獲得了較優(yōu)的設計參數(shù)組合，當魚體剛度約為30MPa、魚體擺動幅值約為0.13l、擺動頻率約為5.16 Hz時具有較好的游動性能，巡游速度達到約1.4l/s。