閆文彥，黃六一，李玉巖，王 剛，鄔騫力，尤鑫星

(中國海洋大學水產(chǎn)學院，青島 266003)

魚類在游動過程中與流體產(chǎn)生非常復雜的相互作用[1]。魚類會借助側線等器官感覺周圍流體的變化，通過身體的自主擺動克服水流作用，以保持合適的魚體姿態(tài)和游泳速度[2-7]。在此過程中魚會消耗大量的能量，從而對其生存和生長產(chǎn)生重要影響[8-9]。研究魚體的水動力學特性，對了解魚類與水域環(huán)境的適應性、個體行為和集群特性、養(yǎng)殖技術和裝備開發(fā)以及仿生學具有重要的意義。

目前國內(nèi)外學者對于魚體自身的肌肉運動、推進機理、游泳能力等方面開展較為豐富研究。Walters等[10]探究了鯖科魚類的身體形態(tài)和游泳性能；周仕杰等[11]對幾種幼魚曲線游泳能力進行了比較研究；Altringham等[12]研究了在模擬游泳條件下肌肉力量的體外輸出模式；Liao等[13]比較了虹鱒在不同直徑圓柱體后方的游泳狀態(tài)和在自由水流中的游泳狀態(tài)，發(fā)現(xiàn)魚類可以從環(huán)境產(chǎn)生的漩渦中獲取能量，以維持其原有的位置；Liao等[14]從能量的角度解釋了采用卡門步態(tài)游泳的鱒魚的神經(jīng)肌肉控制；Lauder等[15-16]分析了魚類游泳的幾種波動推進的經(jīng)典模式；喬曉冬等[17]分析了魚類游動時反卡門渦街的形成和魚體推力的產(chǎn)生機理，結果表明，反卡門渦街的形成是魚體游動時產(chǎn)生推力的根本原因；周澍欣等[18]針對鯉魚的三維模型進行了數(shù)值模擬計算，對其水阻力情況、流場流速分布和壓力分布情況進行分析；Michele等[19]通過模擬肌肉功能來重現(xiàn)魚體游泳的主要特征，比較了不同游泳步態(tài)的性能，得到了頻率、波長、尾振幅等關鍵參數(shù)與所達到的游泳速度或產(chǎn)生的推力之間的關系，從而獲得其游泳性能；田曉潔等[20]針對仿金槍魚三維模型的流場受力進行對比研究，得到仿生機器魚可有效減少其水阻力，從而提高推進效率；王猛等[21]從亞科層面對中國鯉科魚類游泳能力進行了分析和總結；李明昊等[22]構建了仿生虎鯨直行和轉彎運動負載可靠性分析的狀態(tài)函數(shù)，分析了仿生虎鯨的尾鰭設計變量對載荷的影響。之前大多數(shù)研究重點關注的是不同魚類的游泳模式與推進效率關系，為相關仿生學裝備設計提供服務。隨著野生漁業(yè)資源的衰退，人類通過養(yǎng)殖魚類，獲取高質(zhì)量蛋白質(zhì)，已成為漁業(yè)的重要產(chǎn)業(yè)模式。養(yǎng)殖魚類自身的水動力特性與養(yǎng)殖水域的流速、能量消耗和養(yǎng)殖效率以及養(yǎng)殖技術和裝備設計密切相關，而目前針對魚體自身的水動力學特性的研究較少。

虹鱒呈紡錘狀，體略側扁，因其肉質(zhì)蛋白含量高、膽固醇含量極低，具有很高的食用價值、保健價值及可期的市場前景[23-25]，已成為中國水產(chǎn)養(yǎng)殖的重要品種之一，是黃海冷水團深遠海大型鋼結構網(wǎng)箱主要冷水魚養(yǎng)殖品種之一。

本研究以虹鱒為研究對象，依據(jù)其魚體特征參數(shù)構建虹鱒三維模型，使用數(shù)值模擬方法探究虹鱒在不同流速和沖角下的水動力特性、流場流速分布以及魚體表面壓力分布，并通過水槽試驗驗證。研究結果可為虹鱒的個體和群體行為、養(yǎng)殖水域的選擇、養(yǎng)殖技術和裝備開發(fā)以及仿生學研究提供參考。

1 材料與方法

1.1 虹鱒模型

預試驗發(fā)現(xiàn)小體長單個虹鱒模型在水流作用下的受力很小，儀器測量誤差對結果會產(chǎn)生較大影響。本研究選取大規(guī)格虹鱒(全長為84.22 cm，標準體長78.00 cm，體質(zhì)量6.80 kg)，全長和體質(zhì)量參考《中國鮭鱒魚養(yǎng)殖》[26]；體長、體高、頭長、體寬之間的關系參考李思忠等[27]對虹鱒的研究結果進行換算；除尾鰭外，其他部位魚鰭厚度較小，3D打印較難實現(xiàn)，參照相關研究[20]，忽略其他魚鰭。模型基本參數(shù)見表1。

表1 虹鱒魚模型基本參數(shù)Tab.1 Basic parameters of rainbow trout model

使用SolidWorks軟件構建1∶1的虹鱒魚體數(shù)值模型(圖1a)，利用3D打印制作物理模型(圖1b)。3D打印材料選擇1.75 mm的聚乳酸(PLA)，打印精度為層高0.1 mm，壁厚1.0 mm。

圖1 虹鱒模型示意圖Fig.1 Schematic of rainbow trout model

1.2 模型試驗

1.2.1 儀器設備

模型試驗在中國海洋大學循環(huán)動力水槽中進行，水槽觀察段尺寸為4.0 m × 1.0 m × 1.2 m；利用小威龍vectrino點式流速儀測定試驗流速，量程0～2 m/s，精度± 0.5%；穩(wěn)定流速范圍0.2～0.7 m/s；通過六分力儀傳感器測量試驗過程中魚體受到的水動力，量程0～50 kg，精度± 0.3%。

1.2.2 試驗布設

試驗布設如圖2所示，X方向沿水流方向，Y方向與水槽側壁垂直，Z軸垂直向上；通過兩根細螺紋桿(直徑d=10 mm)將虹鱒模型固定在六分力儀上，六分力儀處于觀察段中部，魚體中心線(吻端A尾鰭中心B的連線)距水面28 cm。魚體頭部的朝向-X方向。試驗水溫為20℃。魚體中心線所在的水平面為魚體中心平面，作為后續(xù)研究流場分布的參考面。

圖2 試驗布設圖Fig.2 Experimental layout diagram

1.3 數(shù)值模擬

1.3.1 湍流模型

對于復雜幾何的外部流動的問題，常采用基于雷諾平均N-S(RANS)方程組的k-ε模型。其控制方程為：

(1)

(2)

式中：Ui表示i方向的濾波速度；uj表示j方向的濾波速度；P表示濾波壓力；ρ表示流體密度。

本研究使用Realizablek-ε模型，可以在雷諾應力上保持與真實湍流的一致[28,29]。Realizablek-ε模型優(yōu)化了湍流粘度公式，從渦量擾動量均方根的精確輸運方程推導出ε方程。

Realizablek-ε模型的湍動能及其耗散率輸運方程為:

(3)

(4)

在上述方程中，Gk表示由于平均速度梯度引起的湍動能產(chǎn)生，Gb是由浮力產(chǎn)生的紊流動能；ym可壓速湍流脈動膨脹對總的耗散率的影響。C2和C1ε是常數(shù)，σk，σε分別是湍動能及其耗散率的湍流普朗特數(shù)，在 Fluent中，作為默認值常數(shù)，C1ε=1.44，C2=1.9，σk=1.0，σε=1.2。

1.3.2 計算域和邊界條件

本研究使用的計算軟件為ANSYS 2021版中的Fluent模塊。計算域與試驗水槽觀察段尺度一致，設置為長度L=4.0 m、寬度W=1.2 m、高度H=1.0 m。虹鱒魚模型在計算域中的布設位置與模型試驗的布設相同。

流體選為20℃下等溫不可壓縮水流，利用壓力-速度分離求解器(Semi-Implicit Method for Pressure Linked Equations-Consistent，SIMPLEC)開展穩(wěn)態(tài)計算。重力加速度為9.8 m/s2。計算域左側為速度入口邊界條件，右側為壓力出口邊界條件，其余為對稱壁面邊界條件(圖3)，確保計算域側面邊界不會對魚體水動力分布產(chǎn)生無關干擾。N-S方程(式2)的對流項和擴散項在空間上的離散利用二階迎風格式。收斂標準設置為1×10-4，迭代計算步數(shù)為1 000步。

圖3 計算域及邊界條件Fig.3 Calculation domain and boundary condition

作為數(shù)值計算中的重要環(huán)節(jié)，網(wǎng)格的質(zhì)量直接影響模擬結果的準確性。一般在數(shù)值仿真計算中采用的網(wǎng)格包括結構化網(wǎng)格和非結構化網(wǎng)格，非結構化網(wǎng)格具有網(wǎng)格生成速度快、質(zhì)量好、數(shù)據(jù)結構簡單等優(yōu)點。本研究采用的非結構化網(wǎng)格以六面體單元為主，利用Poly-Hexcore網(wǎng)格劃分方法(圖4)，在魚體邊緣使用多面體網(wǎng)格與六面體網(wǎng)格實現(xiàn)共節(jié)點連接，以達到提升求解效率與精度的目的[30-33]?？紤]到計算精度和計算資源的平衡，本研究劃分了5種網(wǎng)格進行網(wǎng)格無關性檢驗[34-36]。圖5表示0°沖角下0.4 m/s時魚體阻力系數(shù)隨著網(wǎng)格數(shù)的變化。由圖5可知，當網(wǎng)格數(shù)量達到314萬后，隨著網(wǎng)格數(shù)量的增加，阻力系數(shù)基本保持不變。綜合考慮目前計算機能力和計算精度要求，計算中保證網(wǎng)格數(shù)量大于314萬，根據(jù)不同的計算工況，實際使用的網(wǎng)格數(shù)量均處于400～500萬之間。

圖4 網(wǎng)格示意圖Fig.4 Overview of the grids

圖5 網(wǎng)格無關性驗證Fig.5 Mesh independence verification

1.4 工況設置

根據(jù)虹鱒的游泳速度、黃海冷水團常規(guī)流速以及水槽穩(wěn)定流速范圍，選取試驗流速為0.2～0.7 m/s，以0.1 m/s為間隔。沖角設定為來流方向與魚體中心線所在的縱向中剖面的水平夾角。將5條體長約19 cm的虹鱒放在水槽中試驗，從水槽頂部拍攝的圖像，獲取虹鱒在水流中投影線與來流方向的最大夾角為41.8°(圖6)。數(shù)值模擬時的沖角范圍據(jù)此設置為0°～45°范圍，以5°為間隔。

圖6 虹鱒游動過程中魚體與水流的沖角Fig.6 The attack angle between the fish body andthe current during the swimming of therainbow trout

在水槽中保持魚體沖角15°，測定不同來流速度下的魚體所受水阻力，與數(shù)值模擬結果進行對比驗證。在數(shù)值模擬驗證時，為方便對照，將水槽試驗固定模型的2根螺紋細桿一并設置在數(shù)值模型中，以減少誤差。在數(shù)值方法得到驗證基礎上，利用數(shù)值模擬方法開展所有工況的計算。

1.5 數(shù)據(jù)處理方法

魚體的阻力系數(shù)、側向力系數(shù)和雷諾數(shù)計算公式如下：

(5)

(6)

(7)

式中：Fx—水阻力，N；Fy—側向力，N；Cx—阻力系數(shù)；Cy—側向力系數(shù)；Re—雷諾數(shù)；ρ—水介質(zhì)密度，取998.3 kg/m3；S—受流面積(魚體沿體長方向的中心截面)，m2；U—入流速度，m/s；d—特征長度(體長)，m；v—水的運動黏性系數(shù)(20℃水的運動黏性系數(shù)為1.006×10-6m2/s)。

對監(jiān)測點速度無量綱處理，計算公式如下：

(8)

式中：K—速度比；Ui—速度場中任意點的速度；U—來流速度。

2 結果分析

2.1數(shù)值模擬與水槽試驗結果驗證

圖7為15°沖角時虹鱒魚體水阻力實驗值與數(shù)值模擬結果的對比。由圖7可以看出，模型試驗和數(shù)值模擬得到的阻力曲線趨勢一致，均呈冪函數(shù)關系。在流速0.2～0.6 m/s時較為接近，流速0.4 m/s時相對誤差為2.56%，0.7 m/s時，相對誤差為14.02%。由此表明，本文采用的數(shù)值模擬方法具有較高的計算精度，能夠滿足后續(xù)的研究要求。

圖7 15°沖角時虹鱒魚體水阻力實驗值 與數(shù)值模擬結果的對比Fig.7 Comparison of experimental and numerical resultson drag of rainbow trout at 15°

2.2 魚體水動力特性

2.2.1 水動力系數(shù)

魚體的阻力系數(shù)、側向力系數(shù)與雷諾數(shù)關系如圖8和圖9所示。不同沖角條件下，阻力系數(shù)Cx和側向力系數(shù)Cy隨著雷諾數(shù)增大基本保持不變。表明雷諾數(shù)在1.6×105～5.8×105區(qū)間時，同一沖角條件下，水動力系數(shù)處于“自動模型區(qū)”。

圖8 阻力系數(shù)Cx與雷諾數(shù)Re的變化關系Fig.8 Relationship between drag coefficient and Reynolds number

圖9 側向力系數(shù)Cy與雷諾數(shù)Re的變化關系Fig.9 Relationship between side force coefficientand Reynolds number

由于在試驗范圍內(nèi)，在同一沖角條件下，阻力系數(shù)和側向力系數(shù)基本不變，采用阻力系數(shù)和側向力系數(shù)平均值得到圖10為阻力系數(shù)和側向力系數(shù)與沖角關系曲線。由圖可知，在0°～45°沖角范圍內(nèi)，阻力系數(shù)Cx隨著沖角(α)增大呈現(xiàn)指數(shù)增長趨勢,其關系式為Cx=0.025 2 e0.073 1α(R2=0.994 6)；側向力系數(shù)Cy隨著沖角(α)的增大而呈現(xiàn)線性增長，其關系式為Cy=0.015α(R2=0.977 8)。

圖10 水動力系數(shù)與沖角的關系Fig.10 Relationship between hydrodynamic coefficientand inclinations

2.2.2 魚體周圍速度場特性

本研究取吻端與尾叉最深點連線所在的水平面為參考面(見圖2，魚體中心平面)，探究該平面上的速度分布。將魚體周圍流速大于來流速度(無量綱流速大于1)的區(qū)域定義為高流速區(qū)，低于來流流速(無量綱流速小于1)的區(qū)域定義為低流速區(qū)。圖11顯示了20°沖角時不同流速下魚體中心平面周圍無量綱流速分布。

圖11 20°沖角時不同流速下魚體中心平面周圍流場圖Fig.11 Velocity fields around the fish center plane at different velocities at 20°

由圖可知，當沖角固定時，不同流速下高流速區(qū)和低流速區(qū)分布趨勢基本一致，因此，下文以0.6 m/s為例說明由于沖角的變化而引起的速度場變化情況。

不同沖角、流速0.6 m/s時，魚體中心平面周圍流場分布如圖12所示。由圖12可知，在0°沖角時，魚體左右兩側流速對稱分布，在魚體最寬處流速最大，形成兩側的高流速區(qū)；頭部吻端和尾部附近流場流速較小，在頭部前端形成圓形的低流速區(qū)，沿尾柄形成長條形低流速區(qū)。在沖角0°～15°時，迎流面高流速區(qū)出現(xiàn)于魚體最寬處到魚尾之間，說明沖角和魚體結構曲率對表面的流速分布具有協(xié)同影響。隨著沖角的逐漸增大，魚體背流面的高流速區(qū)逐漸向頭部前端移動，迎流面高流速區(qū)呈現(xiàn)逐漸縮小的趨勢，在0°沖角時，迎流面高流速區(qū)占魚體中心平面面積的2.30%，當沖角達到15°時，迎流面高流速區(qū)面積比例下降至1.14%。當沖角大于15°時，背流面高流速區(qū)形態(tài)在魚體頭部穩(wěn)定后，魚體尾鰭兩側開始逐漸形成影響范圍較小的高流速區(qū)，隨著沖角的增大而增大，并沿著背流面逐漸從尾鰭尖端向軀干移動。這可能是由于邊界層在魚體頭部的背流面附近產(chǎn)生分離之后，魚體最寬處到魚尾之間曲率的變化導致了局部高流速區(qū)的產(chǎn)生。當沖角達30°時，迎流面高流速區(qū)消失。在沖角0°～45°時，迎流面低流速區(qū)隨著沖角增大而逐漸擴大，并由魚體前部吻端位置向后擴散，尾部迎流面低流速區(qū)隨著沖角增大而前移，最終低流速區(qū)覆蓋整個迎流面。周澍欣等[18]對鯉魚的三維模型進行了數(shù)值模擬探究，其流場結果與本研究一致，符合紡錘形魚體周圍流場的實際情況。

圖12 不同沖角下魚體中心平面周圍流場分布圖Fig.12 Velocity fields around the fish center plane with varying inclinations

2.2.3 魚體表面壓力分布

流速0.6 m/s，不同沖角條件下，魚體表面壓力分布如圖13所示。0°沖角時，由于對稱于來流方向，魚體迎流面和背流面的壓力分布一致，負壓集中在體寬最大處。負壓是由于魚體的特定體型造成的結果，這與上節(jié)中流體流經(jīng)魚體體寬最大處時，兩側區(qū)域形成高流速區(qū)的結果相吻合，也與田曉潔等[20]對仿金槍魚的研究結果一致。

圖13 不同沖角條件下魚體表面壓力分布云圖Fig.13 Surface pressure distribution nephogram offish with varying inclinations

魚體迎流面正壓區(qū)在沖角小于20°時，主要集中在頭部和尾部，當沖角大于20°時，迎流面正壓區(qū)逐漸貫穿至整個軀干側面。隨著沖角的增大，迎流面上的負壓區(qū)域逐漸減少，負壓區(qū)逐漸向背、腹方向偏移。迎流面上的正壓區(qū)隨著沖角的增大逐漸增大，沖角0°時，正壓區(qū)占迎流面的23.90%，當沖角達到45°時，迎流面的71.66%均為正壓區(qū)。在沖角0°～20°時，魚體最寬處到魚尾之間由于其曲率的變化，壓力幅值逐漸減小，而上文中迎流面高流速區(qū)同樣出現(xiàn)于魚體最寬處到魚尾之間，這說明魚體表面的速度分布和壓力分布息息相關，符合伯努利原理。

魚體背流面正壓區(qū)在沖角小于10°時，只存在于魚體頭部和尾柄部位，正壓區(qū)在魚體背流面的占比大于16.89%，當沖角大于10°時，頭部正壓區(qū)消失，尾柄部位正壓區(qū)隨著沖角的增大而減小，但在任意沖角條件下均為正壓，這也可以通過對照尾端速度的分布得到印證，當沖角達到45°時，背流面上正壓區(qū)僅有1.35%。魚體背流面負壓區(qū)在0°到15°沖角范圍內(nèi)，隨著沖角的增加從最大體周處開始向頭部擴散，當沖角大于15°時，魚從魚體頭部開始，以魚體標準體長(L)為基準，提取魚體中心平面上的體表不同位置(0.1 L、0.25 L、0.5 L、0.75 L和0.9 L)的壓力數(shù)據(jù)(圖14和圖15)，分析迎流面和背流面上魚體表面的壓力變化規(guī)律。由圖可知，隨著沖角的增大，迎流面表面壓力整體上均呈上升趨勢；從魚體位置角度分析，0.25 L之前(魚體頭部)以及0.9 L位置(尾柄)，在所有沖角條件下，表面壓力整體上處于正壓力作用，其中，0.1 L處表面壓力均為最大，說明魚體頭部受到的表面壓力最大，最大壓力出現(xiàn)在沖角30°時，最大值約為194.22 Pa；從沖角的角度分析，沖角25°～45°，各個部位壓力均為正壓，除尾柄位置外，各位置壓力基本上隨著沖角的增加而增加，與圖13中的壓力分布云圖相一致。背流面表面壓力整體上隨著沖角的增加呈線性下降趨勢；沖角30°之后，各個部位表面壓力均為負壓，其中，頭部的壓力變化最大，由0°沖角的正壓102.31 Pa，下降至負壓-345.91 Pa。所有沖角條件下，0.9 L位置(尾柄)的壓力在0 Pa附近，變化不大，與圖13的壓力云圖相吻合。

圖14 迎流面魚體不同部位的壓力變化Fig.14 The pressure in different parts of the fish bodychanges with varying inclinations in stoss side

圖15 背流面魚體不同部位的壓力變化Fig.15 The pressure in different parts of the fish bodychanges with varying inclinations in slip slope

3 結論

本研究依據(jù)虹鱒魚體特征參數(shù)，構建虹鱒魚三維模型，在數(shù)值模擬和水槽試驗驗證的基礎上，利用數(shù)值模擬方法探究了虹鱒個體在不同流況下的水動力特性、魚體中心平面周圍速度場分布以及魚體表面壓力分布情況。在試驗流速(0.2～0.7 m/s)范圍內(nèi)，處于同一沖角條件下時，魚體水動力系數(shù)均處于自動模型區(qū)范圍內(nèi)；在0°～45°沖角范圍內(nèi)阻力系數(shù)Cx隨著沖角增大從0.027增大至0.599，呈現(xiàn)指數(shù)增長趨勢，側向力系數(shù)Cy隨著沖角的增大而呈現(xiàn)線性增長，最大達到了0.71；由于虹鱒魚體本身的體型特征，隨著沖角的增大，速度場中高流速區(qū)和低流速區(qū)呈現(xiàn)一定的移動規(guī)律，隨著沖角的增大，魚體迎流面高流速區(qū)呈現(xiàn)逐漸后移并縮小的趨勢，背流面的高流速區(qū)逐漸向頭部前端移動，背流面高流速區(qū)形態(tài)在魚體頭部穩(wěn)定后，魚體尾鰭兩側開始逐漸形成影響范圍較小的高流速區(qū)，隨著沖角的增大而增大，并沿著背流面逐漸從尾鰭尖端向軀干移動，說明沖角和魚體結構曲率對魚體周圍的流速分布具有協(xié)同影響；魚體結構與表面的壓力分布密切相關，對于魚體的迎流面而言，表面壓力在不同位置上整體呈上升趨勢，其中魚體頭部受壓最大(最大值可達194.22 Pa)；對于背流面而言，背流面表面壓力整體上隨著沖角的增加呈線性下降趨勢，其中，頭部的壓力變化最大(102.31～-345.91 Pa)。以沖角的變化替代魚體的擺動，忽略魚鰭對魚體水動力學特性可能帶來的影響，此外魚體表面不同位置壓力分布和變化規(guī)律與魚體形狀密切相關，需要在今后開展更為深入的研究。

□