水下航行器阻力計算及結構設計

2023-11-13 16:10:38張宇新李鵬魏博秦洪德

應用科技 2023年5期

關鍵詞：模型

張宇新，李鵬,2，魏博，秦洪德

1. 哈爾濱工程大學船舶工程學院，黑龍江哈爾濱 150001

2. 哈爾濱工程大學煙臺研究院，山東煙臺 264006

21 世紀，海洋劃界爭端、海底油氣資源爭端、漁業資源爭端、深海礦產資源爭端層出不窮[1]。隨著人們對海洋探索的愈加深入，海洋裝備也在不斷地改進升級。自主水下機器人（autonomous under-water vehicle，AUV）是一種能夠自主在水下進行海洋科學研究、水下搜索和救援等任務的無人機器人[2]。智能水下機器人的出現，加快了人類探索海洋的速度，由于其經濟費用低、作業效率高、可在復雜的水域環境中持續作業等優點，針對水下機器人各項技術的研究己經成為目前船舶方向的科研工作者們的重點研究內容之一[3]。

其中，AUV 的型線優化一直是海洋工程領域的重要研究方向之一，它很大程度上影響著航行器的性能以及能耗和生產等問題[4]。AUV 型線設計的手段在于通過減小水流的阻力來提高作業性能，從而使AUV 能夠快速準確地完成各種水下任務。除此之外，舵板是一個不可忽視的部件，AUV 的穩定性和快速性可以憑借對舵板的操控達到預期的運動，對于指導AUV 的型線優化、提高其操控性能具有重要意義。對于結構的優化效果可以通過阻力分析來確定，優化的目的之一便是減小水阻。當前，對AUV 的阻力性能進行研究的方法主要包括試驗和計算模擬。通過試驗可以獲得較為準確的AUV 阻力數據，但對試驗設備的要求較高，存在財力與人力上的限制。相比之下，計算流體力學（computational fluid dynamics，CFD）方法具有成本低、可重復性強等優點，是目前應用最廣泛的AUV 阻力計算方法之一[5]。通過建立AUV 的幾何和流場模型，可以模擬AUV 在水中的運動狀態，計算AUV 受到的阻力及其分布情況。在過去的幾十年里，許多學者針對AUV 結構設計優化及數值算法開展了大量的研究。Alvarez 等[6]對水下航行器的最佳船體形狀進行優化研究，采用模擬退火算法來搜索定義航行器形狀的參數設置，以最小化波浪阻力。金碧霞等[7]對某流線型AUV 進行改進，建立了阻力系數計算模型，優化AUV 艏部結構為分體式，增強結構整流作用的同時又能減小水流阻力。胡克等[8]對幾種不同型線回轉體進行了CFD 數值仿真模擬，給出了一些殼體型線優化的建議。戴鵬[9]針對水下航行器的總體設計進行研究，主要研究對象是螺旋槳效率及阻力等性能參數，并用參數優化法優化了航行器的外形。Gao 等[10]提出了一種通過計算流體動力學方法估算流體動力系數的省時方法。Hong 等[11]通過CFD 方法對便攜式AUV 的水動力特性進行研究，建立了AUV 動力學模型，對水動力系數進行估計。

本文采用數值模擬方法，基于計算流體力學知識原理，使用STAR-CCM+軟件對水下航行器模型的整體阻力進行計算，分析航行器表面壓力和周圍流場特性，并對艏部型線以及舵板截面進行減阻優化，確保所選擇的部件為減阻性能更好的一種。

1 數值計算方法

1.1 控制方程與湍流模型

流體運動需要遵循物理守恒定律，這些守恒定律可通過控制方程來進行數學表達。由于本文涉及的流體為絕熱不可壓縮的牛頓流體，沒有能量的交換，所以本文的控制方程為質量守恒方程與能量守恒方程。

質量守恒方程的表達式為

式中： ρ為密度，t為時間，u、v、w分別為速度矢量U在x、y、z等3 個方向上的分量。

對于牛頓流體在x、y、z共3 個方向上的動量守恒方程（也稱Navier-Stokes 方程）為

式中： μ為動力黏度，Su、Sv和Sw為動量守恒方程的廣義源項，p為流體微元上的壓力。

本文中涉及的流體流動屬于湍流的范疇。一般來說，本領域內認為，無論湍流的運動有多么復雜，非穩態的連續方程以及Navier-Stokes 方程對于湍流的瞬時運動仍然是適用的。

本文選擇采用目前工程研究上應用最廣泛的雷諾平均（reynolds average navier-stokes，RANS）方程方法解決湍流問題，RANS 方法首先將滿足動力學方程的湍流瞬時運動分解成平均運動與湍流運動2 部分，然后將脈動運動部分對平均運動的貢獻通過雷諾應力項加以模擬，也就是通過湍流模型來封閉N-S 方程，使之可以被求解得到結果。RANS 方法忽略了密度脈動帶來的影響，但它同時考慮了平均密度的變化。雷諾時均Navier-Stokes 方程如下：

式中i和j取值為1、2、3。

湍動能k和湍動耗散率ε定義為

在標準k-ε模型中，與之相對應的輸運方程為

式中：Gk為湍動能k的第1 產生項，它是由平均速度梯度引起的；Gb為湍動能k的第2 產生項，它是由浮力引起的；YM為在可壓湍流中脈動擴張的貢獻項；G1ε、G2ε、G3ε為一般經驗常數； σk和 σε均為Prandtl 數，它們分別與湍動能k和耗散率 ε對應；Sk和Sε為源項，由用戶定義。

1.2 數值離散方法

本文采用有限體積法對上述控制方程進行離散，離散過程如下：

用 ?表示各物理量，對微分方程在控制體積V內進行積分，即

離散之后可以得到

式中： ?f為 ?在f面上的對流值，(div?)n為(div?)在f面上的法線方向數值，ρfuf?f·Af為通過f面的質量通量，Af為f面的面積向量，Nface為控制體周圍的單元面數量， Γ?為 ?的擴散系數，S?為源項。

離散之后，使用求解壓力耦合方程組的半隱式壓力修正方法，即壓力耦合方程組的半隱式方法（semi-implicit method for pressure linked equations，SIMPLE）對離散后的方程進行求解。

2 阻力性能計算

2.1 計算模型

本文數值模擬所用的水下航行器計算模型如圖1 所示，模型全艇體長8.3 m，艇體最大半徑0.5 m，坐標系采用直角坐標系，坐標原點選為模型頭部頂端，x軸與艇體的中心對稱軸重合，方向與無攻角時來流方向一致，z軸選為豎直向上。

圖1 水下航行器計算模型

根據計算模型確定計算域范圍，AUV 的運動往往決定了計算域的選擇，把重要的流體區域計算進去。數值計算中，計算域的范圍選擇為前端從AUV 艏部向前延伸1.5 倍艇長，后端從AUV尾部向后延伸2.5 倍艇長，左右兩側沿寬度方向左右各延伸2 倍艇寬，上下兩側深度方向各延伸2 倍艇高，形成一個長方體計算區域，如圖2 所示，并通過體積控制進行了網格加密，確保計算精度。

圖2 航行器模型計算域

2.2 邊界條件

計算域劃分之后進行邊界條件的設置，計算區域的入口處設置為速度入口條件（velocity inlet），在區域–邊界中設置相應的速度值，方向沿x軸的負方向；計算區域的出口處設置為壓力出口條件（pressure outlet），此處是完全發展的流動，通過區域內部外推可以得到出流面的流動情況，且應用此方法不會對上游流動產生影響；其他的控制域邊界條件均設置為速度入口條件，但速度大小設置為0，可以模擬AUV 在水下的真實情況。AUV 的表面設置為無滑移的壁面條件（wall），因為計算域設置的足夠大，因此其四周不受AUV艇身的影響；另外，整個計算域設置為流體。

2.3 網格設置

進行流場求解計算前，要將計算區域離散化，即劃分網格。高質量的網格對于進行可靠且準確的計算流體動力學分析是至關重要的，因為網格質量的好壞不僅決定了工況的計算速度，還會對計算結果的準確可靠性產生很大的影響。

在劃分網格的過程中注意到雖然多面體網格可以最大限度保留計算區域的幾何特征，但是其網格質量不高、尺寸太大，在計算中容易帶來誤差。除此之外，切割體網格的收斂速度要快于多面體網格，所以選用切割體網格進行計算。為更好地捕捉航行器附近的流場細節，在航行器圍殼附近生成體控制，以加密此處的網格，網格劃分結果如圖3 所示。

圖3 航行器模型切割體網格

2.4 網格無關性驗證

為確保阻力計算結果的準確性，在航速選擇為7 kn 的情況下，設置了4 組網格數量逐漸增加的工況進行，在除網格數以外的條件均相同的情況下完成各工況阻力計算對網格的無關性進行驗證，驗證結果如表1 所示。

表1 網格無關性

通過計算結果可以得知，隨著網格數量的加密，阻力值逐漸下降，當網格數量為430 萬和242 萬時，網格增加幾乎一倍，但阻力值變化很小，僅減少了0.25%，可以認為結果已經收斂，驗證了網格的無關性。而由于網格的加密，后者的計算時間遠大于前者，因此綜合考慮計算效率與結果的準確性，本文最終選擇的網格數量為2 416 471。

2.5 計算結果及分析

根據上述設置，通過改變邊界的流入速度，分別計算航速為3、5、7、5、9 和11 kn 時，AUV 所受到的航行阻力，計算結果如表2 所示。為清楚地觀察阻力變化趨勢，繪制受力曲線如圖4 所示。其中Fd為阻力，U為航速。

表2 阻力預報數值

圖4 阻力曲線

根據計算結果可知，AUV 的阻力隨著航速的增加逐漸增大，增大趨勢大致呈二次曲線形式，這與文獻[13-14]中計算水下航行器阻力所得到的結果（列于圖5（a）和圖5（b）中）呈現相同的趨勢，驗證了結果的準確性。

圖5 阻力結果趨勢參考

不同航速下水下航行器的流場速度如圖6 所示。由圖6 可知，隨著航行器速度的增大，周圍流場的速度隨之增大，由于物體的存在，航行器首端和尾端的流場速度均有不同程度的減小，并且隨著航速的增大，尾流場速度受影響的區域范圍逐漸增大。在尾流場速度受影響的區域內，從靠近航行器尾端到遠離航行器的方向流場速度從0 逐漸增大，直至與航速一致。

圖6 航行器周圍流場速度

不同航速下水下航行器的表面壓力如圖7 所示。由圖7 可知，航行器的表面壓力隨著航速的增大逐漸增大，并且在航行器艏部出現壓力最大值的情況。在艏部結構設計中應加強強度，以提高航行器艏部的耐壓性。

圖7 航行器表面壓力

3 艏尾部結構設計

文中的水下航行器模型是經過多重結構設計最終確定下來的阻力最小的最優模型，在結構設計過程中進行了艏部型線設計以及舵板剖面確定等內容。

3.1 艏部型線設計

經大量研究，學者們發現了減阻性能較好的水滴型、MYing 型、半橢型以及魚雷型等型線，可用于水下航行器的艏部型線[15]。本文在艏部型線的選擇上共有4 種方案，除了最終選擇的上述改進型艏部型線以外，還提供了水滴型、半橢型、MYing 型3 種艏部型線供選擇，3 種型線的示意如圖8~10 所示。

圖9 半橢型艇艏

圖10 MYing 型艇艏

水滴型艏部數學表達式為

式中：yx為曲線各點處的半徑；xs為軸向位置；D為最大剖面直徑，即平行中段直徑；Ls為艏部的長度；ns為艏部形狀指數，其值的大小表示艏部曲線的豐滿程度，本設計中選擇ns=2.4。

半橢型艏部數學表達式為

MYing 型艏部數學表達式為

式中n為頭部形狀指數，本設計中選擇n=2。

將4 種方案分別計算5 組速度–阻力值并進行對比，確定艏部型線最優結果，計算以及對比結果如表3 所示。為清晰對比4 種艏部型線的阻力結果，將計算結果繪制曲線圖，如圖11 所示。

圖11 不同艇首阻力曲線對比圖

由計算結果可以發現，當航速很低時，4 種型線的阻力結果相近，但隨著航速的增加，半橢型艏部型線的AUV 阻力值上升是最快的，水滴型和MYing 型艏部型線與其相差很小，只略小于半橢型，整體來看，3 種常規艏部型線的阻力結果幾乎一致。而經過特殊優化的艏部型線的阻力結果上升趨勢明顯小于其他3 種，有較好的減阻性能，是更合理的艏部型線選擇，也就是文中第2 節中計算AUV 阻力是選擇的艏部型線。

3.2 舵板設計

除了艏部型線外，舵的選擇也會對AUV 阻力造成影響，由于作業需要以及結構要求，文中水下航行器的舵的結構以及形式已經確定，現在對舵的截面進行確定，使得提高AUV 升力的同時盡可能減少阻力。本節選擇了2 種不同截面的舵板進行升力和阻力性能比較，選擇結果相對較好的一種供航行器使用。2 種舵板分別是按照NACA0012和NACA0020 共2 種翼型作為橫剖面形狀設計的，它們的計算模型如圖12~13 所示。其中，后掠角均選擇20°，舵高選擇0.3 m，舵寬選擇0.45 m。