側噴射流的仿生魚表面減阻優化研究

劉國勇 ，宋中鑫，朱冬梅

（北京科技大學 機械工程學院，北京 10083）

覆蓋著地球71%面積的海洋蘊藏著豐富的魚類、礦產、能源等資源［1］，對于海洋的開發和爭奪成了很多發達國家的戰略重點.仿魚形水下機器人是一種集機械、智能控制于一體的高科技設備，在民用、軍事等科學研究領域均體現出了廣闊的應用前景和巨大的潛在價值［2］.國內外學者越來越重視仿生機器魚的研究與開發，并取得了很多階段性成果.

在基礎理論方面，Gray［3］認為海豚的肌肉所能提供的功率只相當于與它身體相似的剛體模型以15～20節的時速前進時所需功率的1/7，他從能量守恒的角度提出了一個結論：海豚的游動效率遠遠超出了100%.人們將這個結論推廣到整個魚類，稱為“Gray 悖論”.直到今天，這個悖論還激勵著廣大科學工作者以精確的科學方式證明其對錯.Lighthill［4］將空氣動力學的二維機翼理論運用于尾鰭推進的研究，提出了用于分析鲹科魚類推進模式的“細長體理論”，后進一步又提出了“大擺幅細長體理論”［5］.Wu［6］通過研究提出了“二維薄板理論”，將魚看成一個具有彈性的薄板，進一步又提出了“非定常二維波動板理論”［7］，該理論分析了扁平月牙尾魚的推進及優化方式.在應用研究方面，Triantafyllou 等［8］開發并研制出了名為RoboTuna 的機器魚，它是世界上第一臺真正意義上的仿生機器魚，能通過擺動尾柄和尾鰭像魚一樣在水中游動，最高游動速度可以達到 2 m/s.Wang等［9］設計開發了各種類型的機器魚，并將剛體動力學與Lighthill的大幅度細長體理論相結合，提出了由尾鰭驅動的機器人魚的完整動態模型.國內對仿生魚的研究起步較晚，中國科學院自動化研究所研發的多連桿仿生海豚，可實現偏航與俯仰的運動，且首次在水下機器人上實現了混合雜技特技［10］.張藝博等［11］以鱷魚為原型設計了一種可實現水陸兩棲運動的仿生機器人，彌補了國內外的水下機器人機動性能差、復雜地形適應性弱等缺陷.

在減阻方面，仿生減阻、微氣泡減阻、超空泡減阻等各種減阻技術相繼被提出并應用于實踐中，取得了很好的經濟效益，具有很好的科研和應用價值［12］.仿生射流減阻技術是一個新興的研究領域，國內外研究者們對側向射流和逆向射流技術進行了一系列相關的研究［13-14］，發現均有很好的減阻效果.趙剛等［15］基于魚類鰓裂部位仿生射流表面理論分析，對仿生射流表面回轉體進行射流試驗，發現仿生射流表面具有較好的減阻效果，最大減阻率達到10.8%.李芳等［16-17］對流體在多孔仿生射流表面上的流動特性進行了數值模擬，發現減阻效果隨著流速比的增大而增強，最大減阻率為59.02%，并解釋了射流表面減小摩擦阻力的原因及對近壁區邊界層的控制行為.

在仿生魚領域，現有的研究主要針對仿生魚的控制、驅動、力學特性等方面進行研究，對表面減阻的研究相對比較少.此外，當前關于水下表面減阻的研究，主要是基于平板模型或圓筒模型開展，對于流場環境復雜的曲面模型研究較少.本文依據金槍魚的外形建立仿生魚模型，將模型與射流減阻技術相結合，研究射流孔的形狀、高度、位置、高寬比對仿生魚表面減阻的影響，并采用響應面法分析不同因素相互作用時對減阻性能的影響，得出最優減阻方案，對仿生減阻有一定參考意義.

1 模型建立

1.1 模型基本參數

1.1.1 仿生魚模型基本參數

金槍魚適合高速游動，可長時間以30～50 km/h（8.3～13.9 m/s）的速度巡游，瞬間速度最高可達150 km/h，在仿生領域備受關注.本文以金槍魚為仿生對象建立光滑模型，巡游速度范圍取8～14 m/s，使用解析式流線型模仿金槍魚的外形輪廓，忽略金槍魚的魚鰭，仿生魚的幾何參數如圖1 所示.其中L為魚體長度，LE為進流段長度，LR為去流段長度，D為最大截面直徑，DT為尾端截面直徑.本文主要研究射流孔的結構參數對仿生魚模型的減阻影響，定義仿生魚模型的尺寸參數為：L=100 mm，LE=30 mm，LR=70 mm，D=30 mm，DT=3 mm.

圖1 仿生魚的幾何參數示意圖Fig.1 Schematic diagram of the geometric parameters of the bionic fish

1.1.2 射流孔基本參數

設計了前月牙形、后月牙形、前三角形、后三角形4種不同形狀的射流孔，如圖2所示.其中s為射流孔寬度，h為射流孔高度.射流孔來源于對鯊魚鰓部的仿生模擬，通過觀察鯊魚腮部大小，其高度介于形體最大直徑的1/6 至1/3，本文研究的射流孔高度為6～10 mm，高寬比ks=h/s，范圍定為4～10.

圖2 4種形狀的射流孔Fig.2 Four types of jet holes

將射流孔垂直于來流速度方向布置在模型上，如圖3 所示.射流孔到魚首的距離定義為射流孔位置x，本文研究的射流孔位置為5～25 mm.

圖3 射流孔的布置示意圖Fig.3 Schematic diagram of the jet hole arrangement

射流表面減阻率η的計算式為：

式中：FS為光滑表面所受黏性阻力；FJ為射流表面所受黏性阻力.定義相對射流流量Vs：高8 mm 的后三角形射流孔，在射流速度為1 m/s 時，單位時間內的射流流量為V*，相同時間內在任意射流速度下的射流流量為V，則V/V*為該工況下的相對射流流量Vs.從而可以進一步定義單位體積射流減阻率ηs：

式中：ηs為射流表面減阻率η與相對射流流量Vs的比值，表示單位體積流量的射流對減阻率的貢獻.在射流孔大小、射流孔形狀或射流速度不同時，以ηs來比較不同工況下的射流減阻的效率，作為評價射流減阻性能的一個重要指標.

1.2 網格劃分

本文研究的是仿生魚繞流問題，設置大小為Φ400 mm×400 mm 的計算域，為減小工作量，建立 1/2 三維模型.在仿生魚距離魚首5 mm 處，添加高度h=8 mm，高寬比ks=4，射流孔角度為0°（射流孔垂直魚體軸線布置）的前月牙形射流孔.由于模型結構較為復雜，難以對其進行結構網格劃分，對模型進行Tetra/Mixed 非結構網格劃分，并使用Robus（tOctree）劃分策略.這種方法在網格劃分過程中首先生成模型體網格，在這之后生成模型面網格，對復雜結構模型的網格劃分適應能力更強，比較適合用于本節的仿生魚射流減阻模型.

數值模擬計算結果的準確性與網格的數量和密度有著很大關系，而在關注流體表面壓差阻力和摩擦阻力時，需要根據邊界層的實際情況設置合適的邊界層網格高度.本文研究的仿生魚繞流問題為湍流問題，在定義邊界層網格高度時應使第一層網格的y+＜10，以精確求解邊界層流場的物理問題，本文取y+=2，在來流速度為14 m/s 時計算得出近壁面第一層網格的高度設置約為0.003 mm，網格高度增長比率設置為1.001，全局網格尺寸設置為10 mm，仿生魚模型表面的最大網格尺寸設置為0.2 mm，網格總數為400 萬，如圖4 所示，下文模型的網格劃分在此基礎上進行.其中，網格劃分使用ICEM 15.0，仿真計算在Fluent 15.0中進行.

圖4 射流模型的網格劃分Fig.4 Grid division of jet model

Quality 是CFD 常用的網格質量評估標準之一，對于不同類型的網格其采用不同的衡量方式，對于三角形或四面體網格，其代表著計算高度與每一條邊的長度比值.其表示網格的整體質量，數值為0～1，數值為1時最好，通常認為0.3以上時網格質量較好.圖5 為三維模型網格質量，從圖5 中可以看出，網格的Quality大部分處于0.38以上，網格質量較好，滿足計算要求.

圖5 三維模型網格質量Fig.5 3D model mesh quality

1.3 邊界條件設置

本文研究主要關注仿生魚表面阻力，其中摩擦阻力（黏性阻力）與邊界層處流體的剪切應力有關，需要精確求解邊界層中的黏性底層，因此湍流模型選擇SSTκ-ω的模型，使用低雷諾數方程求解.求解器選擇基于壓力的分離求解器，分離算法選擇SLMPLEC 算法，離散格式選擇二階迎風離散格式（Second Order Upwind），收斂精度設置為十萬分之一，流場介質設置為液態水、不可壓縮流體，密度為998.2 kg/m3，黏度為0.001 003 kg（/m·s），通過計算可得出仿生魚流場周圍雷諾數為8.0e+5～1.4e+6.流場入口為速度入口（Velocity-inlet），出口為壓力出口（Pressure-outlet），并將壓力出口的表壓（Gauge Pressure）設置為默認為0.流場外邊界、仿生魚表面以及旋轉壁面設置為無滑移壁面邊界條件（Wall）.流場入口和出口的湍流黏度比（Turbulent Viscosity Ratio）設置為默認的5%，其余設置保持默認.1/2 模型計算域如圖6 所示.下文的模型分析在此基礎上進行.

圖6 1/2模型計算域Fig.6 Calculation domain of one-half model

1.4 網格無關性驗證

仿真計算受到計算機計算能力的限制，因此需要網格獨立性檢驗，來獲得合適的網格尺寸和數量，以保證仿真模擬計算的精度和效率.采用3 種方案對仿生魚模型進行結構網格劃分，網格數量分別為200 萬、400 萬、800 萬.保持其他設置相同，分別將網格導入Fluent 中進行計算.3 種方案下計算得出的阻力系數如表1 所示，其中阻力系數變化比為網格數量200 萬和400 萬時計算得到的結果分別與網格數量為800萬時的計算結果相比較得到的.

表1 3種網格數量下的阻力系數計算結果及對比Tab.1 The calculation results and comparison of resistance coefficient under three kinds of grid numbers

由表1 中可以看出，網格數量從400 萬增加到800萬時，各項阻力變化比在1%以內，表明網格數量為400萬的量級能夠滿足計算精度要求.

2 主流場速度及射流孔結構參數對減阻的影響

2.1 主流場速度對減阻的影響

在仿生魚距離魚首5 mm處，添加射流孔高度h=8 mm，形狀為前月牙形的射流孔.在主流場速度（來流速度）v=8～14 m/s（間隔1m/s 取值），射流速度VJ=0.2～1.2 m/s（間隔0.2 m/s 取值）的情況下，分析主流場速度v對減阻效果的影響.

圖7 為不同主流場速度下總阻力減阻率與速度比VJ/v的關系.從圖中可以看出，在射流速度VJ與主流場速度v的比值保持不變的情況下，不同主流場速度下的減阻率比較接近.這說明，合理調整射流速度VJ的大小，在主流場速度較大時也能夠獲得較高的減阻率.

圖7 不同主流場速度下減阻率與速度比的關系Fig.7 Relationship between drag reduction rate and velocity ratio under different main flow field velocities

圖8 為不同主流場速度下仿生魚模型的表面單位體積射流減阻率.從圖中可以看出，在主流場速度相同時，射流速度VJ越大，仿生魚模型的表面單位體積射流減阻率越高.在相同射流速度VJ下，主流場速度v越大，單位體積射流減阻率越低.

圖8 不同主流場速度下的單位體積射流減阻率Fig.8 Drag reduction rate of unit volume jet under different main flow field velocities

2.2 射流孔結構參數對減阻的影響

在主流場速度v=14 m/s，射流速度VJ=0.2～1.2 m/s（間隔0.2 m/s 取值）的情況下，分別單獨改變射流孔形狀、射流孔高度、射流孔位置、射流孔高寬比這4個參數，分析其對減阻效果的影響，仿真分析結果如圖9所示.

由圖9（a）可以看出，隨著射流速度VJ的增大，不同形狀射流孔的單位體積射流減阻率呈現出先增大后趨于穩定的趨勢，在相同射流速度下，后三角形的射流孔減阻效果更好.由圖9（b）可以看出，在同一高度下隨著射流速度的增大，單位體積射流減阻率逐漸增大，隨著射流孔高度h的增加，不同射流速度下的單位體積射流減阻率呈現出先增大后減小的趨勢.由于單位體積射流減阻率表示單位體積流量的射流對減阻率的貢獻，改變射流孔位置并不會改變單位體積流量，因此選取總阻力減阻率來評價不同射流孔位置的減阻效果，如圖9（c）所示.從圖中可以看出，在射流孔位置x相同的情況下，減阻率隨著射流速度VJ的增加而增大，且射流孔越靠后，總阻力減阻率隨射流速度增加而增大的幅度就越小.在射流速度相同的情況下，射流孔越靠前，總阻力減阻率就越大，且射流速度越大，減阻效果越明顯.由圖9（d）可以看出，當射流孔高寬比ks=4 時，仿生魚模型的表面單位體積射流減阻率隨射流速度的增大而增加，并趨于穩定值；當ks＞4 時，仿生魚模型的表面單位體積射流減阻率隨射流速度的增大而減小最終趨于穩定.在射流速度VJ相同的情況下，高寬比ks為6 時單位體積射流減阻率最大.

3 參數優化設計

在傳統的設計過程中，往往需要考慮許多參數，而每個參數又需要足夠的變量來提供樣本數據，這樣工作量是隨著參數數量的增多指數型增長的，因此往往采用優化設計的技術進行參數優化.本文擬采用響應面法分析射流孔各個結構參數在交互作用時對減阻的影響，并求出最優設計方案.

3.1 實驗設計

影響仿生魚表面減阻的結構參數主要有射流孔形狀、射流孔高度、射流孔位置及射流孔高寬比等.根據Box-Behnken 中心組合試驗設計原理，本文使用Design-Expert 軟件進行三因素三水平試驗，射流孔高度、射流孔位置、射流孔高寬比為自變量，仿生魚模型的總阻力為響應值，其試驗水平及仿生魚的總阻力分別見表2 和表3.確定相關因素不變量為：射流孔形狀為后三角形，主流場速度v=14 m/s，射流速度VJ=1.2 m/s.

表2 響應面分析與因素水平表Tab.2 Response surface analysis and factor level table

表3 響應面實驗設計及結果Tab.3 Response surface experimental design and results

3.2 回歸方程建立與分析

將表3 中的數據導入Design Expert 軟件中進行多元回歸擬合，建立回歸方程如下：

式中：R表示模型總阻力的預測值.

對回歸方程進行顯著性分析，結果如表4所示.

表4 模型總阻力回歸方程分析Tab.4 Analysis of the model total resistance regression equation

模型的F值為50.24，P值為0.000 2，表明該模型是顯著的，可信度較高.回歸方程中的h、x、hks、xks、h2、x2較為顯著，其余三項不夠顯著，但整體來看，該二次回歸方程模型比較顯著，與實際值能夠較好的擬合，可以使用此模型對仿生魚模型的總阻力大小進行預測.

3.3 響應面分析

分析比較射流孔高度、射流孔位置、射流孔高寬比3個參數對總阻力的影響，使一參數固定不變，比較其余兩參數對總阻力的影響，結果如圖10～圖12所示.

圖10 射流孔高度與射流孔位置對仿生魚總阻力影響Fig.10 Effect of jet hole height and jet hole position on total drag of bionic fish

等高線的形狀能夠反映交互效應的情況，橢圓形表示交互程度強，圓形表示交互效應弱.

由圖10 可以看出，在射流孔的高寬比固定時，高度h和位置x交互作用等高線圖呈現橢圓形，表明二者的交互作用顯著.由圖10（a）可知，仿生魚的總阻力隨著高度h和位置x的增加，呈現出先上升后下降的趨勢，且位置x的上升幅度大于高度h，說明射流孔位置對總阻力的影響較大.當射流孔高度為 10 mm，射流孔位置為5 mm時，模型預測的仿生魚總阻力達到最小值2.547 82 N.當射流孔高度固定時，位置x與高寬比ks對總阻力的影響如圖11 所示.由圖11（a）可知，隨著位置x和高寬比ks的增大，總阻力呈現出先上升后下降的趨勢，且位置x的上升幅度大于高寬比ks，說明射流孔位置對總阻力的影響大于射流孔高寬比.由圖11（b）可知，射流孔位置x與射流孔高寬比ks的交互作用顯著.當射流孔位置為 5 mm，射流孔高寬比為4 時，模型預測的仿生魚總阻力達到最小值2.496 82 N.如圖12 所示，由圖12（a）可知，隨著高度h和高寬比ks增加，總阻力呈現出先上升后下降的趨勢，但這一趨勢并不明顯，高度h的上升幅度相對而言要大于高寬比ks，說明射流孔高度對總阻力的影響大于射流孔高寬比.當高度為6 mm，高寬比為10 時，模型預測的總阻力達到最小值2.601 06 N.

圖11 射流孔位置與射流孔高寬比對仿生魚總阻力影響Fig.11 Effect of the position of jet hole and the ratio of height and width of jet hole on the total drag of bionic fish

圖12 射流孔高度與射流孔高寬比對仿生魚總阻力影響Fig.12 Effect of jet hole height and jet hole ratio on total drag of bionic fish

利用Design-Expert軟件分析優化，通過系統爬坡模擬尋優法，以總阻力最小為目標得到的最優方案為：射流孔高度為6 mm、射流孔位置為5 mm、射流孔高寬比約為4，此時模型預測可達到的仿生魚總阻力最小值為2.483 46 N.采用優化后的方案，進行仿真驗證，重復3 次，取其平均值得到的仿生魚總阻力為2.510 21 N，對應的減阻率為6.49%.

4 結論

本文以金槍魚為仿生對象，建立了仿生魚的光滑形體模型和添加射流孔的模型，采用數值模擬的方法，分析了主流場速度、射流孔結構參數等因素發生改變時對仿生魚表面阻力產生的影響.使用響應面法分析了不同射流孔結構參數在共同作用時對仿生魚表面減阻的影響，得出了最優減阻方案.

得出的主要結論如下：

1）當來流速度增加時，需要增大射流速度，使得射流速度與來流速度的比值較大，才能獲得較高的減阻率.

2）在相同射流速度下，后三角形的射流孔減阻效果更好；隨著射流孔高度h的增加，不同射流速度下的單位體積射流減阻率呈現出先增大后減小的趨勢；在所給的范圍內，射流孔越靠前，總阻力減阻率就越高，且射流速度越大，減阻效果越明顯；當高寬比ks為4時，單位體積射流減阻率隨射流速度的增大而增加，當ks＞4 時，單位體積射流減阻率隨射流速度的增大而減小最終趨于穩定.

3）在射流孔形狀選用后三角形、射流孔高度為6 mm、射流孔位置為5 mm、高寬比為4 時，能夠達到更好的減阻效果，通過仿真驗證得到了6.49%的減阻速率.