高鵬騁 劉冠杉 黃橋高 ,2) 潘 光 馬云龍
* (西北工業大學航海學院,西安 710072)
? (中國船舶集團有限公司系統工程研究院,北京 100094)
魚類集群運動的現象隨處可見,集群游動能幫助魚群中各個體有效規避捕食者,提升捕食效率,增加繁殖機會等[1-3].此外,許多學者認為魚類可以通過高度組織化的群體運動來保存能量,獲得水動力優勢,提高個體的運動性能.
大多數研究認為,魚群游動的水動力優勢是由于渦流假說和槽道效應[4].其中,渦流假說認為,魚群中后魚置身于前魚的尾渦中,通過傳入的尾渦結構的誘導速度降低了魚和水流間的相對速度,后魚同時還能借助前魚尾渦提高游動效率.槽道效應指出,通過增加相鄰魚之間在游泳方向的水流速度來降低魚和水流間的相對速度,同時魚類還可以借助相鄰魚渦街誘導出的流動來提高游動效率.
近年來,隨著計算方法以及數值仿真軟件的發展,學者們對魚群水動力性能進行了大量仿真計算,并得到了許多有意義的結論.Deng 等[5]采用改進的浸入邊界法對菱形魚群中的3 條魚單元進行了數值模擬.結果表明,位于前一列兩條魚之間的后魚可以利用上游魚脫落的反卡門渦街來提高推進效率以及降低功耗.隨后,Chung[6]在文獻[5]的基礎上進行了深入研究,給定了3 條魚最優的位置及運動參數設置,即兩條前魚橫向間距為0.4 倍弦長,進行反相波動運動,后魚的前緣距前魚后緣0.2 倍弦長.Tian 等[7]通過有限元法數值求解不可壓縮的Navier-Stokes 方程來研究子母魚對的游泳性能和渦結構,計算結果表明母魚和小魚在串聯和交錯排列中都受益,通過改變它們的相位差和相速度,可增加推力和降低功耗.目前,還有不少學者利用強化學習算法來對先導-跟隨群游結構中的跟隨魚進行運動控制,其中先導魚的運動方式固定,跟隨魚通過強化學習來調整運動方式,以實現跟隨魚利用先導魚的渦流來提升推進效率[8-9].
王亮[10]對二維仿真魚群游動進行了數值模擬,指出當魚群中各個體存在體形差異時主要利用側向水動力(即槽道效應)來提高推進效率;當個體體形相近時,主要利用前魚尾渦提高效率.Li 等[11]對串聯、并聯、菱形和矩形4 種編隊的魚群游動進行了數值模擬.結果表明,效率受到尾流和側向壓力的影響,其中尾流主要影響推力,橫向功率損失受到側向壓力影響.Lin 等[12-15]對二維鰻式游動魚群開展了系列研究,探究了相位差、間距對串聯魚群及傾斜布置魚群的水動力性能影響.
隨著電機驅動技術,傳感器技術以及流場顯示技術的發展,一些學者開始設計搭建實驗裝置來研究魚群游動的節能機理[16-21].Dewey 等[16-17]與Boschitsch 等[18]分別實驗研究了兩個柔性撲翼在并聯與串聯布置時的水動力性能,并利用DPIV 系統進行流場記錄及顯示.實驗結果表明,對于并聯結構,兩個翼形的推進特性相同: 同相拍動時,效率提升而推力減小(相對于單個翼形);反相拍動時,推力增加而效率基本不變.對于串聯結構,前翼的推力和效率在分離距離較大時與單個翼形一致,而后翼由于受到前翼尾流的影響,其動力學特性與分離距離和相位差相關.
當前國內外關于集群水動力特性研究廣泛集中于二維翼型而非三維生物的研究.并且,研究主要集中在二維平面的群游,而沒有考慮到集群中各單體沿垂向分布的情況,然而通過Vicsek 等[22]的研究表明生物在實際游動過程中往往采取上下錯落分布隊形.此外,目前的研究對象大多為尾鰭擺動推進魚群,忽略了對采用滑撲結合推進的鰩科生物的研究.
蝠鲼屬鰩科生物,有著扁平形的身體、較大的展弦比,巡游速度在0.25~0.47 m/s[23].蝠鲼采用撲動與滑翔相結合的推進方式[24],效率可高達89%,高于鰻科、鲹科等水中生物[25-27],在機動性、隱蔽性等方面具有超凡的優越性.為彌補現有研究的不足以及探究雙蝠鲼垂向間距和攻角(α)對各個蝠鲼水動力性能的影響,本文針對沿垂向分布的雙蝠鲼進行滑翔運動時的水動力特性開展研究,其中垂向間距選取了0.25TL(TL為體厚),0.5TL,0.75TL,1TL4 種工況,在每種工況下進行了攻角變換,變化范圍為?8°~8°(每2°取一個特征攻角),以期為仿蝠鲼滑撲一體航行器在滑翔時的隊形設置提供參考數據.
采用逆向工程技術建立蝠鲼數值計算模型,通過測量實際物體的尺寸,獲得物體的三維點數據,再通過點數據構建三維曲線,進一步構建三維曲面,進而重構實物的CAD 模型[28],如圖1 所示.其中展長(SL)為2900 mm,體長(BL)為1800 mm,最大體厚為350 mm.

圖1 蝠鲼模型Fig.1 Manta ray model
我們將4 種間距排布統一展示在圖2 中,計算域尺寸根據實際排布方式進行適應性調整,其中保證入口與蝠鲼間距4 倍體長,出口與蝠鲼間距7 倍體長,展向壁面距蝠鲼3 倍展長,上、下壁面分別距上下蝠鲼3 倍體厚.邊界條件設置如下: 入口邊界設為速度入口,出口邊界設為壓力出口,壁面設為自由滑移壁面,蝠鲼表面設為靜止無滑移壁面.

圖2 計算域設置Fig.2 Calculation domain setting
本文在4 組垂向間距(0.25TL~1TL)以及9 組攻角(?8°~8°)組合下進行了仿真計算,探究了集群滑翔時,垂向間距和攻角對集群中各單體和集群系統的水動力特性的影響.
根據文獻[29],雙體蝠鲼滑翔狀態下的系統平均阻力(Ddouble)按照下式計算
其中,Dup和Ddown分別代表了上、下蝠鲼所受到的阻力.
將蝠鲼集群滑翔過程中的受力無量綱化
其中,D,L為蝠鲼所受的阻力、升力,ρ為流體密度,U為來流速度,BL為蝠鲼體長.
本文借助ICEM 進行結構網格繪制,由疏到密繪制了4 套網格(網格數量分別為1.5 × 106,3 × 106,4.5 × 106,6 × 106)進行單體仿真計算,以期確定出在保證計算精度前提下兼顧計算速度的網格節點設置(如圖3(a)所示).圖3(b)展示了在4 套網格下,蝠鲼單體在0°攻角下以0.5 m/s 速度滑翔的阻力系數(CD)、升力系數(CL)數值仿真結果.可以看出當單體網格設置為4.5 × 106時,計算結果精確度已經達到最優,為節省計算成本,提高計算速度,在后續的雙體蝠鲼仿真計算中我們沿用了此時的節點設置.

圖3 網格細節及無關性驗證Fig.3 Grid detail and independence verification
本文借助 FLUENT 軟件進行計算,湍流模型設置為SSTk-ω模型,采用 SIMPLEC 算法對連續方程中的壓力和速度的耦合,收斂殘差設置為1.0 × 10?5.
采用此方法對三維漸變NACA0009 翼型在雷諾數Re=1.0 × 106時不同攻角下的阻力系數(CD)、升力系數(CL)進行了仿真計算,并與文獻[30]的實驗數據進行比對,如圖4 所示.由圖4 可知: 數值仿真所得到的CD和CL與實驗值接近,確保了數值計算方法的可靠性.

圖4 方法驗證Fig.4 Method validation
本節展示了雙蝠鲼以0.5 m/s 速度滑翔時,不同垂向間距下阻力(升力)系數隨攻角的變化情況.研究發現,間距和攻角對各阻力系數產生較大影響,但對升力系數,尤其是系統平均升力系數影響較小.因此,本節按照垂向間距劃分小節,在每小節中結合壓力云圖對各阻力系數進行分析,在垂向間距為0.25TL時結合蝠鲼表面壓力云圖探究升力系數變化.
本小節對垂向間距0.25TL的雙蝠鲼集群滑翔進行了數值模擬,集群中各單體攻角保持同步變化,變化范圍為?8°~8°.圖5 展示了上(下)方蝠鲼阻力(升力)系數、雙蝠鲼系統平均阻力(升力)系數以及單體蝠鲼阻力(升力)系數隨攻角變化的改變情況.由圖5(a)可以看出,系統平均阻力均大于單體滑翔時所受到的阻力.但當雙蝠鲼以負攻角滑翔時,下方蝠鲼受到阻力遠低于單體滑翔,尤其是在?8°攻角時,幾乎不受到阻力;當以正攻角滑翔時,上方蝠鲼受到的阻力減小,下方蝠鲼受到阻力增加.與此同時,我們還發現當攻角為?2°時,上下兩蝠鲼受到的阻力幾乎相同.由圖5(b)可以看出,當滑翔攻角小于?2°時,系統平均升力大于單體滑翔時的升力;當滑翔攻角大于?2°時,系統平均升力小于單體滑翔時的升力;當滑翔攻角為?2°時,系統平均升力與單體滑翔升力接近.此外,我們還發現上方蝠鲼所受升力始終小于下方蝠鲼.

圖5 0.25TL 垂向間距時不同攻角下的阻力/升力系數Fig.5 Drag/lift coefficient at different attack angles at 0.25TL vertical distance
為進一步探究阻力(升力)出現變化的原因,我們將壓力/速度云圖繪制在了圖6 中,其中壓力取值范圍為?300~150 Pa,速度取值范圍為0~0.8 m/s,后續所有云圖中取值范圍保持一致.從圖6(a)和圖6(b)可以看出,無論以何種攻角進行集群滑翔時,蝠鲼外圍壓力場分布相似,即在蝠鲼的頭部及尾部附近各存在著一個高壓區,高壓區的面積及壓力大小幾乎一致.當蝠鲼集群以不同攻角集群滑翔時,集群中各單體出現阻力差異的原因在于各單體的表面壓力出現變化.由圖6(a)可以看出,下方蝠鲼頭部的中下部分處于低壓區,而上方蝠鲼頭部完全處于高壓區,這必然導致下方蝠鲼的阻力減小,上方蝠鲼的阻力增加.由圖6(b)可以看出,當以正攻角進行集群滑翔時,上方蝠鲼頭部的中上部分處于低壓區,下方蝠鲼頭部完全處于高壓區,因此下方蝠鲼的阻力增加,上方蝠鲼的阻力減小.接下來,結合圖6(c)解釋上下兩蝠鲼升力表現,從圖6(c)中可以看出,上下兩蝠鲼上表面壓力分布基本相似,下方蝠鲼頭部壓力略高于上方蝠鲼;下方蝠鲼下表面壓力明顯高于上方蝠鲼,尤其是在蝠鲼的中腹部,上方蝠鲼存在著明顯的低壓區.又因為升力受表面壓力影響較大,因此下方蝠鲼所受升力總是大于上方蝠鲼.由圖6(d)可以看出兩蝠鲼中間區域速度最高,尾部速度最低,頭部存在小面積低速區,因此在蝠鲼尾部及頭部出現高壓區,在兩蝠鲼中腹部周圍出現低壓區,這和壓力云圖相契合.通過結果統計,我們發現雙蝠鲼沿垂向分布集群滑翔時各升力系數幾乎不受間距和攻角的影響,在后續小節中為避免重復,僅結合壓力云圖對阻力表現進行分析.與此同時,我們發現隨著間距的增加,集群滑翔對速度場影響逐步減弱,我們后續僅在間距為1TL時展示速度云圖.

圖6 0.25TL 垂向間距時不同攻角下的壓力/速度云圖Fig.6 Pressure/velocity diagram at different attack angles at 0.25TL vertical distance
同3.1 節,圖7 展示了當垂向間距為0.5TL時,上(下)方蝠鲼阻力(升力)系數、系統平均阻力(升力)系數以及單體蝠鲼阻力(升力)系數隨攻角變化的改變情況.由圖7(a)可以看出,系統平均阻力仍大于單體滑翔時所受到的阻力,但相較于0.25TL間距時,有所下降.上(下) 方蝠鲼的阻力變化趨勢同0.25TL一致.由圖7(b)可以看出,系統平均升力變化不大,上下兩蝠鲼升力差值減小.
同樣,我們將?8°,4°時的壓力云圖繪制在了圖8 中,壓力云圖分布情況同間距為0.25TL時大致相同,從圖8(a)可以看出,在以負攻角滑翔時,上方蝠鲼頭部的中下部開始出現低壓區,這將導致上方蝠鲼的阻力出現下降,這與圖5 和圖7 的阻力系數曲線圖相符.當以正攻角滑翔時,上下蝠鲼的表面壓力相較于間距為0.25TL時沒有明顯變化,因此在推力曲線上也幾乎一致.

圖7 0.5TL 垂向間距時不同攻角下的阻力/升力系數Fig.7 Drag/lift coefficient at different attack angles at 0.5TL vertical distance

圖8 0.5TL 垂向間距時不同攻角下的壓力云圖Fig.8 Pressure diagram at different attack angles at 0.5TL vertical distance
同前,圖9 展示了當垂向間距為0.75TL時,上(下) 方蝠鲼阻力(升力) 系數、系統平均阻力(升力)系數以及單體蝠鲼阻力(升力)系數隨攻角變化的改變情況.由圖9(a)可以看出,系統平均阻力更加接近單體滑翔時所受到的阻力,上、下方蝠鲼的阻力差值也減小,說明垂向間距的增大,集群中單體間相互影響降低.同樣,我們將?8°,4°時的壓力云圖繪制在了圖10 中,從圖10(a)可以看出,在以負攻角滑翔時,上方蝠鲼頭部的中下部低壓區面積擴大,使得上方蝠鲼阻力進一步降低(見圖9).從圖10(b)可以發現以正攻角滑翔時,下方蝠鲼頭部開始出現低壓區,使得下方蝠鲼的阻力減小(見圖9).由此可以看出,隨著垂向距離的增加,集群中各單體間的相互作用減弱,系統平均阻力和單體滑翔阻力逐步接近.由圖9(b)可以看出,系統平均升力變化不大,上下兩蝠鲼升力差值進一步減小.

圖9 0.75TL 垂向間距時不同攻角下的阻力/升力系數Fig.9 Drag/lift coefficient at different attack angles at 0.75TL vertical distance

圖10 0.75TL 垂向間距時不同攻角下的壓力云圖Fig.10 Pressure diagram at different attack angles at 0.75TL vertical distance
同前,圖11 展示了當垂向間距為1TL時,上(下)方蝠鲼阻力系數、系統平均阻力系數以及單體蝠鲼阻力系數隨攻角變化的改變情況.由圖11(a)可以看出,當攻角在?8°,?4°時,系統平均阻力與單體滑翔時所受到的阻力幾乎保持一致,上、下方蝠鲼的阻力差值也進一步減小,集群中單體間相互影響進一步削弱.由圖11(b)可以看出,系統平均升力變化不大,隨著垂向距離的增加,上下兩蝠鲼的升力差值也隨之減小.我們將?8°,4°時的壓力/速度云圖繪制在了圖12 中,從圖12(a)和圖12(b)可以看出,當以負攻角滑翔時,上方蝠鲼頭部的中下部低壓區面積擴大,使得上方蝠鲼阻力進一步降低,類似的,當以正攻角滑翔時,下方蝠鲼頭部的低壓區面積也有所增加,因此此時下方蝠鲼的阻力也隨之減小.由圖12(c)可以看出隨著間距的增大,上下兩蝠鲼的速度云圖幾乎保持一致,集群滑翔對集群中各單體的速度場影響減弱.


圖11 1TL 垂向間距時不同攻角下的阻力/升力系數Fig.11 Drag/lift coefficient at different attack angles at 1TL vertical distance

圖12 1TL 垂向間距時不同攻角下的壓力/速度云圖Fig.12 Pressure/velocity diagram at different attack angles at 1TL vertical distance
本文借助Fluent 軟件對沿垂向分布的雙蝠鲼集群滑翔狀態進行了數值仿真,得到了升力、阻力系數曲線,并結合壓力云圖進行分析,得到了如下結論.
(1)系統平均升力受垂向間距影響較小,當滑翔攻角小于?2°時,系統平均升力大于單體滑翔時的升力;當滑翔攻角大于?2°時,系統平均升力小于單體滑翔時的升力;當滑翔攻角為?2°時,系統平均升力與單體滑翔升力接近.下方蝠鲼升力總大于上方蝠鲼,但隨著垂向間距增加,上下兩蝠鲼升力差距減小.
(2)雙蝠鲼沿垂向分布在攻角范圍為?8°~8°進行集群滑翔時未能發現系統平均阻力低于單體滑翔時所受阻力的情況,但隨著垂向間距的增加,單體間相互影響降低,系統平均阻力趨近于單體滑翔阻力.
(3)系統平均阻力雖未能發現減阻效果,但集群中單體能從集群滑翔中獲得減阻收益.當雙蝠鲼以負攻角集群滑翔時,下方蝠鲼獲得減阻收益,且垂向間距越小,獲得的減阻效果越明顯;當以正攻角集群滑翔時,則是上方蝠鲼獲得減阻收益.
(4) 雙蝠鲼集群滑翔過程中各單體阻力(升力)表現的差異性主要來源于蝠鲼頭部、中腹部壓力區的變化,垂向間距的變化改變了壓力區分布,從現有結果可以看出隨著垂向間距的進一步加大,上(下)蝠鲼壓力分布會趨于一致.