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

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

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

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

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

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

1 數(shù)值計算方法

1.1 控制方程與湍流模型

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

質(zhì)量守恒方程的表達式為

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

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

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

本文中涉及的流體流動屬于湍流的范疇。一般來說，本領(lǐng)域內(nèi)認為，無論湍流的運動有多么復雜，非穩(wěn)態(tài)的連續(xù)方程以及Navier-Stokes 方程對于湍流的瞬時運動仍然是適用的。

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

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

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

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

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

1.2 數(shù)值離散方法

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

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

離散之后可以得到

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

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

2 阻力性能計算

2.1 計算模型

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

圖1 水下航行器計算模型

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

圖2 航行器模型計算域

2.2 邊界條件

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

2.3 網(wǎng)格設置

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

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

圖3 航行器模型切割體網(wǎng)格

2.4 網(wǎng)格無關(guān)性驗證

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

表1 網(wǎng)格無關(guān)性

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

2.5 計算結(jié)果及分析

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

表2 阻力預報數(shù)值

圖4 阻力曲線

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

圖5 阻力結(jié)果趨勢參考

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

圖6 航行器周圍流場速度

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

圖7 航行器表面壓力

3 艏尾部結(jié)構(gòu)設計

文中的水下航行器模型是經(jīng)過多重結(jié)構(gòu)設計最終確定下來的阻力最小的最優(yōu)模型，在結(jié)構(gòu)設計過程中進行了艏部型線設計以及舵板剖面確定等內(nèi)容。

3.1 艏部型線設計

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

圖9 半橢型艇艏

圖10 MYing 型艇艏

水滴型艏部數(shù)學表達式為

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

半橢型艏部數(shù)學表達式為

MYing 型艏部數(shù)學表達式為

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

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

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

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

3.2 舵板設計

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

圖13 NACA0020 型舵板

NACA0012 幾何表達式為

式中：t=0.12，b為舵寬（取0.45），x、y為橫縱坐標。

NACA0020 幾何表達式為

式中：t=0.20，b為舵寬（取0.45），x、y為橫縱坐標。

2 種舵板的阻力和升力的計算結(jié)果如表4 所示。為更清晰化比較，將結(jié)果作圖對比，如圖14和圖15 所示，其中，U表示航行速度，F(xiàn)d表示阻力，L表示升力。

表4 2 種舵板不同航速下計算結(jié)果

圖14 2 種舵板阻力對比

圖15 2 種舵板升力對比

由計算結(jié)果可知，隨著航速的提高，2 種舵板的阻力和升力都隨之增加，但NACA0012 型舵板的阻力增加更緩慢，并且升力增加更迅速，是更好的舵板截面選擇。因此最終本水下航行器選擇安裝的是NACA0012 型舵板。

4 結(jié)論

通過STAR-CCM+軟件對本文水下航行器進行結(jié)構(gòu)優(yōu)化設計以及阻力計算，結(jié)果顯示：

1）航行器整體阻力值整體隨著航速的增大而增大，大致呈二次曲線形式。不同的艏部型線對航行器的阻力有較大的影響，使用優(yōu)化型艏部型線可以較其他型線使阻力增漲的更緩慢。因此本文選擇的優(yōu)化型艏部型線較其他常規(guī)性艏部型線有更好的減阻效果。

2）航行器周圍流場速度隨著航速的增大而增大。航行器使其首尾端的流場速度有不同程度的減小，且隨著航速的增大，尾流場速度受影響的區(qū)域范圍逐漸增大。

3）航行器表面壓力同樣隨著航速的增大而增大，并且在航行器艏部達到最大值，也就是說航行器的最艏部是整體結(jié)構(gòu)中受力最大的，在結(jié)構(gòu)設計時應著重考慮。

4）對于舵板的選擇，NACA0012 型舵板不僅比NACA0020 型舵板可提供更大的升力，并且阻力值也更小，是水下航行器更優(yōu)的舵板選擇。