龐永杰，王亞興，楊卓懿，高 婷

（哈爾濱工程大學(xué)水下智能機器人技術(shù)重點實驗室，黑龍江哈爾濱150001）

Myring型回轉(zhuǎn)體直航阻力計算及艇型優(yōu)化

龐永杰，王亞興，楊卓懿，高 婷

（哈爾濱工程大學(xué)水下智能機器人技術(shù)重點實驗室，黑龍江哈爾濱150001）

在水下智能機器人（AUV）方案設(shè)計階段，針對如何在主尺度和巡航速度確定的前提下得到直航阻力最小的回轉(zhuǎn)體形狀這個問題，以Myring型回轉(zhuǎn)體為研究對象，利用計算流體力學(xué)（CFD）方法計算其在水中做勻速直航運動時的阻力并與試驗值比較以驗證準(zhǔn)確性。分別對二維和三維2種網(wǎng)格形式，標(biāo)準(zhǔn)和加強2種壁面函數(shù)的計算精度和計算效率進行了對比，并在此基礎(chǔ)上探討單獨改變艏部或艉部形狀時阻力的變化規(guī)律。再利用EXCEL、ICEM、FLUENT建立優(yōu)化平臺，基于多島遺傳算法在限制條件范圍內(nèi)對Myring型回轉(zhuǎn)體參數(shù)進行全局尋優(yōu)，尋找不同流速時阻力最小的Myring型回轉(zhuǎn)體，為AUV設(shè)計階段的阻力性能預(yù)報和艇型優(yōu)化提供參考。

水下智能機器人；回轉(zhuǎn)體；阻力；模型試驗；多島遺傳算法；流體力學(xué)；網(wǎng)格生成；壁面函數(shù)

水下智能機器人（autonomous underwater vehicles，AUV）作為完成多種水下任務(wù)的重要工具，無論在軍事還是科研領(lǐng)域都有廣闊的前景。足夠的續(xù)航力和巡航速度是AUV完成預(yù)定任務(wù)的前提。目前AUV不能達到預(yù)定續(xù)航力和速度的原因主要是低估了AUV的航行阻力或AUV本身直航阻力過大。本文通過經(jīng)驗公式，計算流體力學(xué)方法（CFD）對AUV航行中所受的阻力進行預(yù)估并參考實驗數(shù)據(jù)進行驗證，尋找滿足預(yù)定設(shè)計要求且能夠有效降低直航阻力的Myring型回轉(zhuǎn)體解決方案。

1 概述

1.1 Myring型回轉(zhuǎn)體

回轉(zhuǎn)體AUV的艇型一般由公式給出艏艉形狀曲線，并根據(jù)需要決定是否使用平行中段，常用的艏艉形狀有 Myring型［1］、Nystrom 型［2］、卡克斯型［2］、水滴型［3］等，邱遼原［4］在潛艇粘性流場數(shù)值模擬研究中也提出了一種艇型。這些艇型都可以通過改變公式中的參數(shù)值來獲得不同的形狀。所有回轉(zhuǎn)體形式中Myring型使用得較多，著名AUV如Remus，MAYA都使用此形式，本文對Myring型回轉(zhuǎn)體進行系統(tǒng)分析，尋找以Myring方程為基礎(chǔ)的回轉(zhuǎn)體直航阻力的確定方法，為AUV設(shè)計中預(yù)估艇體阻力和有效降低總阻力提供參考。Myring給出的艏部形狀方程為

艉部形狀方程為

式中：a為艏部長度，b為平行中體長度，c為艉部長度，d為中體直徑，x是長軸上點到艏部頂點的距離，r為x點處的半徑，n和θ分別是控制艏艉曲線飽和程度的參數(shù)。n和θ越大艏艉越飽滿。圖1給出了a＝25，b＝15，c＝30，d＝10的時候不同n和θ值下艏艉的形狀。

圖1 不同n值和θ值時Myring方程給出的艏艉形狀Fig.1 Structure diagram of inductive energy storage shape of revolution body when different n and θ are adopted

1.2 回轉(zhuǎn)體阻力計算

獲得回轉(zhuǎn)體直航阻力的方法主要有經(jīng)驗公式、CFD方法和模型試驗。Barros將計算導(dǎo)彈和魚雷航行阻力的DATCOM方法借鑒到AUV艇體阻力計算中［5］，Prestero在計算Remus艇體阻力［6］時引用了Triantafyllou的成果。經(jīng)驗公式最大的優(yōu)點就是計算迅速，但是計算精度不夠，只能在設(shè)計初期為設(shè)計者提供一定的參考。隨著計算機技術(shù)的飛速發(fā)展，CFD方法計算AUV水動力性能也得到大范圍的應(yīng)用，受計算機水平的影響，直接數(shù)值模擬（DNS）和大渦模擬（LES）尚處于發(fā)展階段，利用時均化的方法求解Navier-Stokes方程（RANS）已經(jīng)廣泛應(yīng)用于計算流體力學(xué)領(lǐng)域，并發(fā)展出多種湍流模型。CFD方法較經(jīng)驗公式求解更為準(zhǔn)確，耗費的時間也相對長一些，但與模型試驗相比耗時要短很多，且成本也遠(yuǎn)小于后者。模型試驗是獲得AUV水動力參數(shù)最直接有效的方法，但成本也最高，而且受試驗設(shè)備精度和尺度效應(yīng)影響，仍會存在誤差。

2 CFD計算和模型試驗

2.1 CFD計算

CFD計算回轉(zhuǎn)體直航阻力有三維和二維2種計算方法。前者直接模擬真實的三維空間流場環(huán)境，利用體網(wǎng)格求解三維空間下的N-S方程并對回轉(zhuǎn)體表面積分得到阻力值，后者則使用面網(wǎng)格，將回轉(zhuǎn)體的對稱軸作為旋轉(zhuǎn)軸邊界條件，直接求解二維空間上的N-S方程，解得回轉(zhuǎn)體表面沿長度方向上的壓力分布，利用回轉(zhuǎn)體的軸對稱性求解阻力，這種方法的原理可以參考OpenFoam關(guān)于楔形體計算的相關(guān)理論［7］。

二維和三維結(jié)構(gòu)化網(wǎng)格見圖2。通過計算驗證二維和三維模型CFD計算之間的誤差不超過2%。

圖2 二維和三維結(jié)構(gòu)化網(wǎng)格Fig.2 3D and 2D structured mesh

對流場采用有限體積法求解RANS方程，湍流模型選擇在近壁區(qū)算法更穩(wěn)定且精度更好的剪切應(yīng)力輸運k-ω模型［8］，即SST k-ω模型。為使計算更加準(zhǔn)確，收斂更快，使用結(jié)構(gòu)化網(wǎng)格對流場進行劃分，為避免流域邊界對流場產(chǎn)生影響，流場應(yīng)足夠大，整個外域采用圓柱體。根據(jù)Barros的建議，流場的總長度取艇長的15倍，直徑為艇體直徑的25倍［9］。雖然SST k-ω模型在近壁區(qū)的計算上有優(yōu)勢，但仍需借助壁面函數(shù)法將壁面上的物理量與湍流核心區(qū)內(nèi)的物理量聯(lián)系起來。不同的壁面函數(shù)對第一層網(wǎng)格高度的無量綱參數(shù)y＋要求不同，y＋的表達式為

式中：△y是第一層網(wǎng)格高度，μτ是壁面摩擦速度，ρ是水的密度，τω為壁面切應(yīng)力，

標(biāo)準(zhǔn)壁面函數(shù)［10］要求y＋值應(yīng)盡量滿足30 ＜y＋＜60。如果采用加強壁面函數(shù)［11］，則應(yīng)盡量滿足y＋＝1，最大不能超過5。表1給出了選擇不同的y＋時得到的阻力值。此時來流速度為1.5 m／s。試驗值為9.79 N。

表1 y＋值對模型阻力計算的影響Table 1 Effect of y＋to model drag calculation

雖然使用標(biāo)準(zhǔn)壁面函數(shù)較加強壁面函數(shù)誤差偏大，但仍不超過2%。考慮到后續(xù)的優(yōu)化工作需要進行大量的重復(fù)計算，本文仍然選擇標(biāo)準(zhǔn)壁面函數(shù)。

2.2 模型試驗

為獲得不同長度艏艉配合的回轉(zhuǎn)體阻力性能，設(shè)計并制作一組玻璃鋼模型，包括可拆卸的艏部與艉部各兩個以及一個平行中段。模型中段直徑d＝280 mm，長度b＝734 mm。兩組艏部參數(shù)分別為a＝280 mm，n＝2和a＝504，n＝1.8；兩組艉部參數(shù)分別為c＝504，θ＝38°和c＝784，θ＝27°。兩組艏和兩組艉與平行中段分別進行組合，組合形式如圖3，從上到下依次為組合1到組合4。

圖3 模型試驗4種組合形式Fig.3 4 types of model assembly

將4種組合分別在水平型循環(huán)水槽中進行試驗，改變來流速度，利用六分力天平測量各種組合的受力情況，并提取其中的軸向力進行比較。試驗?zāi)Ｐ鸵妶D4。圖5顯示的是組合1在循環(huán)水槽中進行試驗。試驗中測力天平布置在回轉(zhuǎn)體內(nèi)部軸線上，利用兩個劍桿固定測力天平和回轉(zhuǎn)體模型，為將劍桿對回轉(zhuǎn)體周圍流場的影響降到最小，劍桿底部的直徑僅為10 mm，約為模型尺度的5%，而且本試驗中回轉(zhuǎn)體不運動，故劍桿對回轉(zhuǎn)體的軸向上的受力影響很小。

圖4 試驗?zāi)Ｐ虵ig.4 Test model

圖5 組合1在進行試驗Fig.5 Assembly 1 is carrying out experiment

2.3 結(jié)果對比

圖6給出組合1和組合4通過CFD計算和模型試驗EXP獲得的回轉(zhuǎn)體阻力值，組合2和組合3結(jié)果相同。CFD計算與模型試驗獲得D的阻力值基本一致，但在來流速度過高時模型試驗值明顯高于CFD計算值，這主要是因為循環(huán)水槽造流能力的限制，生成流速大于1.5 m／s時伴隨產(chǎn)生過高的湍流，流場不穩(wěn)定使阻力增加。

圖6 組合1和組合4阻力曲線Fig.6 Drag curves of assembly 1 and 4

3 利用CFD計算優(yōu)化艇型

要利用CFD方法進行艇型優(yōu)化需要不斷改變艇型參數(shù)，建立模型并劃分網(wǎng)格，再利用CFD程序迭代求解RANS方程獲得阻力值。三維模型進行一次CFD計算約需要60～80萬體網(wǎng)格，利用4個CPU來迭代2 h以上才可能收斂得到阻力值，二維模型一次CFD計算則僅需9～12萬面網(wǎng)格，用一個CPU經(jīng)過5 min的計算即可得到阻力值。以1 000次計算為例，前者需要80 d，后者則僅需4 d，但兩者計算結(jié)果相差并不大。本文利用二維模型進行艇型優(yōu)化。

3.1 艏部和艉部的優(yōu)化

雖然對于回轉(zhuǎn)體來說中段的長度越小對減小阻力越有利，但AUV的設(shè)計中考慮到總體布置需要，有時會確定平行中段的長度b，再考慮改變艏部和艉部形狀以達到減小直航阻力的目的。即在b和d都固定的前提下討論a和n（或b和θ）變化對直航阻力的影響。本文取試驗參數(shù)作為建模參考，即b＝734 mm，d＝280 mm，由于本試驗的目的是為AUV的實際設(shè)計提供依據(jù)，取來流速度為AUV的預(yù)定巡航速度3 kn（1.5 m／s）。

圖7（a）給出c＝784 mm，θ＝27°固定不變，a＝d、2d、3d、4d時n從0.3增加到4.9過程中CFD計算模型阻力值的變化。從中可以看出在計算范圍內(nèi)不論n取何值，總阻力都隨a值的增大而增加，除a＝d時n值的增加對總阻力的影響不大外，在其他情況下n值的增加都會使總阻力增加，且在n＜2時阻力的增加量特別明顯。這一結(jié)論和Myring的結(jié)論［1］有一定的差別，主要是因為Myring的試驗是在Re＝107的條件下進行的，而本試驗Re＝3×106，此時摩擦阻力在總阻力中的比例較大，而a越大，濕表面積越大，致使摩擦阻力變大。因而Myring的試驗結(jié)論［1］直接運用于航速較低的AUV艇型優(yōu)化具有局限性。

圖7（b）給出a＝504 mm，n＝1.8保持不變，c＝2d、3d、4d、5d時θ從1°增加到35°過程中CFD計算模型阻力的變化。c＝2d時的阻力θ的增加對總阻力幾乎沒有影響，越大的c值的總阻力隨θ變化幅度越大。

圖7 阻力隨a和n，c和θ變化曲線Fig.7 Drag curves of n when different a is adopted，c when different θ is adopted

3.2 固定總?cè)莘e的艇型優(yōu)化

3.2.1 計算模型

AUV設(shè)計中可能根據(jù)任務(wù)需求和設(shè)備大小情況先確定需要的總?cè)莘e，并給出艇體總長度和直徑的比值所存在的范圍。針對AUV直航阻力的艇型優(yōu)化也假定回轉(zhuǎn)體的體積和細(xì)長比已經(jīng)確定，利用Myring方程定義回轉(zhuǎn)體曲線尋找阻力最優(yōu)解，即在V和d／l范圍已知的前提下考查a、b、c、n和θ不停變化時直航阻力的變化。為與試驗數(shù)據(jù)進行對比，并考慮為水下航行器的設(shè)計提供參考，取l＜10d，0.09 m3＜V＜0.15 m3對，Myring型回轉(zhuǎn)體進行優(yōu)化。

利用優(yōu)化程序集成EXCEL、ICEM和FLUENT，集成原理圖見圖8，首先由優(yōu)化組件提供a、b、c、n和θ的初始值，利用計算器作為條件限制組件保證輸出值的正確性，再用內(nèi)置在EXCEL中的VBA程序讀入初始值并輸出型值文件，ICEM執(zhí)行RPL文件利用型值構(gòu)建二維計算模型并輸出非結(jié)構(gòu)網(wǎng)格文件，F(xiàn)LUENT執(zhí)行JOU文件導(dǎo)入網(wǎng)格進行迭代，迭代初步穩(wěn)定后據(jù)AUV的邊界層y＋值和流場的速度梯度對網(wǎng)格進行自適應(yīng)后再計算，自適應(yīng)前后網(wǎng)格的變化見圖9，待迭代穩(wěn)定后提取阻力值。優(yōu)化組件根據(jù)求得的阻力值重新生成a、b、c、n和θ并重復(fù)上述步驟，直到尋找到最優(yōu)解。優(yōu)化過程中使用非結(jié)構(gòu)網(wǎng)格而不是結(jié)構(gòu)網(wǎng)格主要是因為非結(jié)構(gòu)網(wǎng)格能更好的適應(yīng)AUV型線的不斷變化，而且初步計算后自適應(yīng)新增的加密網(wǎng)格本身也是非結(jié)構(gòu)的。本算例的幾何結(jié)構(gòu)非常簡單，所以非結(jié)構(gòu)網(wǎng)格計算結(jié)果與結(jié)構(gòu)網(wǎng)格計算結(jié)果幾乎完全相同。只是在收斂速度方面前者比較差，故在優(yōu)化過程中增加了FLUENT迭代的次數(shù)。優(yōu)化中來流速度分別取v＝0.5、1.0、1.5、2.0 m／s。

圖8 優(yōu)化過程Fig.8 Optimizing process

圖9 自適應(yīng)前后網(wǎng)格Fig.9 Mesh before and after self-adaption

優(yōu)化中改變各參數(shù)但保證a＞200 mm，b＞50 mm，c＞100 mm，d＞100 mm，0.6＜n＜3，5°＜θ＜35°不斷變換參數(shù)值并提取迭代計算獲得的阻力值。為達到尋找最優(yōu)解并避免得到局部最優(yōu)解的目的，優(yōu)化方法選擇全局探索方法中的多島遺傳算法［12］（multi-island and genetic algorithm，MIGA）優(yōu)化二維模型。為保證尋優(yōu)過程的有效性并充分利用計算機資源，本文MIGA的種群規(guī)模取30，進化代數(shù)20，交叉概率取1，變異概率0.01，遷移概率0.1，每次計算用時6 min，在每一個來流速度下共計算2 000次。

3.2.2 優(yōu)化結(jié)果

圖10是來流速度為v＝0.5 m／s時利用多島遺傳算法，以阻力值最小和體積值最大同時作為優(yōu)化目標(biāo)得到的計算結(jié)果，從結(jié)果看在某一體積值下，無論回轉(zhuǎn)體的形狀如何變化，都不會小于某一阻力值Dmin，且 Dmin值隨體積的增加近似呈線性規(guī)律增加。圖 11顯示的是來流速度分別為0.5、1.0、1.5、2.0 m／s時Dmin隨體積增加的變化規(guī)律。雖然速度越大Dmin隨體積變化的越快，但在這幾個來流速度下的阻力最優(yōu)艇型基本上都如圖12所示。平行中段很小甚至沒有，艏部長度和艉部長度的比值a／c約為4，總長度與直徑的比值l／d約為8，n不超過2.5，θ不超過30°。

圖10 v＝0.5 m／s時隨體積增加回轉(zhuǎn)體阻力計算結(jié)果Fig.10 Drag of different volumes when v＝0.5 m／s

圖11 不同速度下最小艇體阻力隨體積變化圖Fig.11 Minimal drag curves at different speeds

圖12 阻力最優(yōu)艇型Fig.12 Best revolution body shape with minimal drag

4 結(jié)束語

本文以Myring型回轉(zhuǎn)體為研究對象，利用CFD方法研究了改變回轉(zhuǎn)體艏或艉的長度和飽滿程度時對回轉(zhuǎn)體直航阻力的影響。搭建了一個優(yōu)化平臺，在固定艇體總?cè)莘e的前提下利用多島遺傳算法尋找改變回轉(zhuǎn)體艏、中、艉3部分在總長度中的比例以及外形參數(shù)時阻力最優(yōu)的方案，同時得到了一種阻力最優(yōu)艇型。對水下航行器的設(shè)計工作具有很好的參考價值。

［1］MYRING D F.A theoretical study of body drag in subcritical axisymmetric flow［J］.Aeronautical Quarterly，1976，27（3）：186-194.

［2］杜月中，閔健，郭字洲.流線型回轉(zhuǎn)體外形設(shè)計綜述與線型擬合［J］.聲學(xué)技術(shù)，2004，23（2）：93-101.

DU Yuezhong，MIN Jian，GUO Zizhou.A review and mathematical formulation of shape design of streamlined bodies of revolution［J］.Technical Acoustics，2004，23（2）：93-101.

［3］吳寶山，潘子英，夏賢，等.水滴型回轉(zhuǎn)體帶尾導(dǎo)管的水動力特性研究［J］.船舶力學(xué)，2003（6）：54-59.

WU Baoshan，PAN Ziying，XIA Xian，et al.Investigation of the hydrodynamic characteristics of body of revolution with stern ring-wing［J］.Journal of Ship Mechanics，2003（6）：54-59.

［4］邱遼原.潛艇粘性流場的數(shù)值模擬及其阻力預(yù)報的方法研究［D］.武漢：華中科技大學(xué)，2006：79-84.

QIU Liaoyuan.Numerical simulation of the viscous flow over the submarine and method research on predicting its viscous resistance［D］.Wuhan：Huazhong University of Science and Technology，2006：79-84.

［5］De BARROS E A，PASCOAL A，De SA E.Investigation of a method for predicting AUV derivatives［J］.Ocean Engineering，2008，35：1627-1636.

［6］PRESTERO T.Verification of a six-degree of freedom simulation model for the REMUS autonomous underwater vehicle［D］.California：University of California，1994：25-27.

［7］Open Foam user guide，version 2.1.1［EB／OL］.（2012-12-05）.http：／／www.openfoam.com.

［8］王福軍.計算流體動力學(xué)分析［M］.北京：清華大學(xué)出版社，2004：202-204.

［9］De BARROS E A，DANTAS J L D.Effect of a propeller duct on AUV maneuverability［J］.Ocean Engineering，2012，42：61-70.

［10］SALIM M S，CHEAH S C.Wall y＋strategy for dealing with wall-bounded turbulent flows［C］／／Proceedings of the International Multiconference of Engineers and Computer Scientists.Hong Kong，2009（2）：2165-2170.

［11］MULVANY N J，CHEN Li，TU Jiyuan，et al.Steady-state evaluation of two equation rans turbulence models for high-Reynolds number hydrodynamic flow simulations［R］.Victoria：DSTO Platform Sciences Laboratory，2004：14-15.

［12］鐘毅芳，陳柏鴻，王周宏.多學(xué)科綜合優(yōu)化設(shè)計原理與方法［M］.武漢：華中科技大學(xué)出版社，2007：64-68.

Direct route drag calculation and shape optimization of Myring shape axisymmetric revolution body

PANG Yongjie，WANG Yaxing，YANG Zhuoyi，GAO Ting
（State Key Laboratory of Science and Technology on Autonomous Underwater Vehicle，Harbin Engineering University，Harbin 150001，China）

When main dimensions and cruising speed have been predefined，how to get an axisymmetric revolution body with minimal drag force at the conceptual design stage of autonomous underwater vehicles（AUVs）has been a problem.This paper analyzed drag force of Myring shaped axisymmetric revolution body moving underwater in a direct route，at constant speed，by CFD method and compared with the experiment value to check the veracity.2D or 3D mesh，standard or enhanced wall treatment function was adopted respectively to figure out the accuracy and efficiency.And on this basis，the change law of drag when changing merely length and shape of the nose or tail was discussed.Based on multi-island and genetic algorithm （MIGA），an optimization platform was constructed with EXCEL，ICEM and FLUENT to search for the best axisymmetric revolution body with minimal drag force at different speed when volume is predefined.This is a good reference for calculating and minimizing total drag of AUVs.

autonomous underwater vehicles（AUV）；revolution body；drag；model test；multi-island and genetic algorithm（MIGA）；fluid dynamics；mesh generation；wall function

10.3969／j.issn.1006-7043.201304028

U674.941

A

1006-7043（2014）09-1093-06

http：／／www.cnki.net／kcms／doi／10.3969／j.issn.1006-7043.201304028.html

2013-04-08. 網(wǎng)絡(luò)出版時間：2014-08-29.

國防科學(xué)工業(yè)委員會基礎(chǔ)研究基金資助項目（9140C270305120 C2701，A2420110001）.

龐永杰（1955-），男，教授，博士生導(dǎo)師.

龐永杰，E-mail：pangyongjie＠hrbeu.edu.cn.