何 隆,張 亞,李世中,王芙蓉,袁澤慧,?；壑?/p>

(中北大學(xué) 機(jī)電工程學(xué)院， 太原 030051)

1 引言

浩瀚的海洋蘊(yùn)藏著各類豐富的資源，隨著開發(fā)利用海洋資源的進(jìn)程加快，對(duì)各類無人航行器的要求日益提升。其中，水下滑翔機(jī)能夠依靠調(diào)整自身重浮力的大小，以“之”字形的運(yùn)動(dòng)軌跡在不同的剖面收集海洋信息，這種航行器能耗較低，但其運(yùn)動(dòng)單一，抗干擾力差，工作效率低，無法勝任形式多樣的探測(cè)任務(wù)[1]。水下多旋翼航行器通過推進(jìn)器差速可以實(shí)現(xiàn)各種復(fù)雜的運(yùn)動(dòng)，適應(yīng)力強(qiáng)，但由螺旋槳旋轉(zhuǎn)產(chǎn)生的噪聲與振動(dòng)無法避免，且受限于能源供應(yīng)而無法長(zhǎng)時(shí)間、遠(yuǎn)距離地進(jìn)行作業(yè)[2]。其次，潛航器的水動(dòng)力性能直接取決于其外形，許多學(xué)者為提高AUV的運(yùn)動(dòng)性能，對(duì)外形的優(yōu)化設(shè)計(jì)作了深入研究，如：吳利紅等[3]通過CFX軟件計(jì)算了不同攻角和漂角工況下滑翔機(jī)所受的位置力，并對(duì)滑翔機(jī)展開了線型優(yōu)化；趙凱凱等[4]利用CFX軟件對(duì)滑翔蛇形機(jī)器人進(jìn)行了水動(dòng)力仿真，求解了水動(dòng)力參數(shù)，并采用矩形平板翼對(duì)翼型進(jìn)行了優(yōu)化；李天博等[5]以最大化滑翔機(jī)的升阻比為優(yōu)化目標(biāo)，建立了Kriging代理模型，采用EGO算法對(duì)聯(lián)翼式滑翔機(jī)外形進(jìn)行了優(yōu)化設(shè)計(jì)；楚燈旺等[6]通過FLUENT軟件分析了4種浮標(biāo)的外形方案，建立了殼體外形阻力模型，采用基于懲罰函數(shù)的遺傳算法優(yōu)化了殼體線型參數(shù)；張代雨等[7]通過CAESES軟件對(duì)翼身融合滑翔機(jī)進(jìn)行全參數(shù)化建模，以升阻比為優(yōu)化目標(biāo)，并基于多種組合優(yōu)化算法提高了外形的水動(dòng)力性能；高啟升等[8]提出了基于組合加權(quán)響應(yīng)面法的多目標(biāo)優(yōu)化方法，以折衷規(guī)劃法對(duì)子目標(biāo)進(jìn)行歸一化處理，對(duì)梯形肋骨耐壓結(jié)構(gòu)進(jìn)行了多目標(biāo)優(yōu)化設(shè)計(jì)。

在之前的研究基礎(chǔ)上，提出了一種四螺旋槳驅(qū)動(dòng)與重浮力調(diào)整裝置相結(jié)合的混合驅(qū)動(dòng)無人潛航器整體結(jié)構(gòu)方案，通過FLUENT軟件分析了AUV在不同攻角和速度下的水動(dòng)力系數(shù)，利用拉丁超立方采樣試驗(yàn)設(shè)計(jì)確定了樣本點(diǎn)，建立了二階響應(yīng)面模型，以升阻比最大和質(zhì)量最小為優(yōu)化目標(biāo)，采用多目標(biāo)遺傳算法對(duì)響應(yīng)面模型進(jìn)行優(yōu)化求解，得到了AUV的優(yōu)化設(shè)計(jì)方案，為AUV水動(dòng)力外形進(jìn)快速優(yōu)化設(shè)計(jì)提供一種高效的方法。

2 結(jié)構(gòu)方案

AUV采用模塊化設(shè)計(jì)，其結(jié)構(gòu)組成可分為殼體、多傳感器陣列單元、浮力調(diào)節(jié)單元、重心調(diào)節(jié)單元、吊裝背鰭、外掛設(shè)備單元、主控與驅(qū)動(dòng)單元、推進(jìn)器單元、尾翼和電力電子單元，結(jié)構(gòu)組成如圖1所示。

圖1 AUV結(jié)構(gòu)組成示意圖Fig.1 AUV structure composition

其中浮力調(diào)節(jié)單元通過步進(jìn)電機(jī)帶動(dòng)活塞移動(dòng)改變單元艙內(nèi)進(jìn)水量實(shí)現(xiàn)凈浮力的調(diào)節(jié)功能。重心調(diào)節(jié)單元通過2個(gè)步進(jìn)電機(jī)改變艙內(nèi)物體的相對(duì)位置實(shí)現(xiàn)重心的遷移。通過上述2個(gè)單元的配合可以實(shí)現(xiàn)AUV低功耗的運(yùn)行模式，按照“之”字軌跡在海洋中緩慢運(yùn)動(dòng)。推進(jìn)器單元由4個(gè)獨(dú)立的推進(jìn)器組成，通過螺旋槳差速可實(shí)現(xiàn)多自由度的運(yùn)動(dòng)。多傳感器陣列單元包含水聲定位模塊、聲納成像模塊、深度測(cè)量傳感器、GPS等組成。電力電子單元包含航姿模塊，無線電信標(biāo)模塊，電源管理模塊等組成，為便于搬運(yùn)，在重心調(diào)節(jié)單元艙外殼處焊接有背鰭。AUV結(jié)構(gòu)的基本參數(shù)如表1所示。

表1 結(jié)構(gòu)的基本參數(shù)

3 水動(dòng)力參數(shù)計(jì)算

AUV在進(jìn)行作業(yè)時(shí)易受到諸如海浪、洋流等各種環(huán)境因素的影響，其水動(dòng)力性能則是判定AUV設(shè)計(jì)好壞的重要依據(jù)，也是實(shí)現(xiàn)精準(zhǔn)控制、精準(zhǔn)作業(yè)的基礎(chǔ)[9]。設(shè)計(jì)時(shí)AUV艏部曲線AB采用雙參數(shù)橢圓型曲線，艉部曲線CD采用圓錐曲線設(shè)計(jì)，整體線型簡(jiǎn)單且加工成本低，如圖2所示。其中回轉(zhuǎn)半徑D為84 mm，中間段長(zhǎng)度L為1 230 mm。

圖2 AUV線型示意圖Fig.2 AUV Linetype

為便于計(jì)算結(jié)果收斂，提升計(jì)算速度，在不改變AUV整體外形結(jié)構(gòu)的前提下省略了對(duì)水動(dòng)力性能影響較小但易造成復(fù)雜網(wǎng)格的零件，通過建模軟件CREO建立了AUV的簡(jiǎn)化模型，并采用右手慣性坐標(biāo)系統(tǒng)，以AUV在靜態(tài)時(shí)的浮心作為坐標(biāo)原點(diǎn)，其中沿原點(diǎn)向艏部為x軸正方向，沿原點(diǎn)向左為y軸正方向，垂直于機(jī)身向上為z軸正方向。AUV的水動(dòng)力模型和坐標(biāo)系如圖3所示。

圖3 水動(dòng)力模型示意圖Fig.3 Hydrodynamic model

3.1 數(shù)值計(jì)算方法驗(yàn)證

為確保數(shù)值計(jì)算結(jié)果的可靠性，以上海交通大學(xué)制作的“海鷗”號(hào)[10]為例對(duì)文中所用到的計(jì)算方法的有效性進(jìn)行驗(yàn)證，“海鷗”號(hào)模型及其網(wǎng)格劃分結(jié)果如圖4所示。求解模型在0°攻角下分別以0.2 m/s、0.3 m/s、0.4 m/s、0.5 m/s的速度直航時(shí)的阻力值，并將其與水池拖曳試驗(yàn)結(jié)果進(jìn)行對(duì)比，如圖5所示。二者數(shù)值相差較小，證明了所運(yùn)用的數(shù)值計(jì)算方法的有效性。

圖4 驗(yàn)證模型與網(wǎng)格劃分結(jié)果示意圖Fig.4 The verified the model and meshing result

3.2 數(shù)值計(jì)算與結(jié)果分析

AUV在水中的航行速度一般不超過4 m/s，屬于低速狀態(tài)，航行中主要受到水的粘性阻力，且隨航行速度的減小其粘性阻力占總阻力的比重將會(huì)增大，因此AUV的水動(dòng)力性能主要體現(xiàn)在水阻力計(jì)算與分析方面。借助FLUENT軟件對(duì)AUV在不同攻角和航速下的阻力系數(shù)、升力系數(shù)展開了計(jì)算。

圖5 數(shù)值驗(yàn)證對(duì)比曲線Fig.5 Comparison curve of numerical verification

為了更好地模擬低湍流雷諾數(shù)流動(dòng)，避免在k-ε湍流模型中重增加復(fù)雜的阻尼函數(shù)，采用RANS方法進(jìn)行數(shù)值求解，使用k-ωSST湍流模型。設(shè)置流體域入口為速度入口，流體域出口為壓力出口，出口壓力為0，AUV表面設(shè)置為無滑移壁面，四周流域邊界設(shè)置為自由滑移壁面。將計(jì)算域分成兩部分，將廣闊的海域簡(jiǎn)化成15L×8L×8L的矩形域，將AUV所在位置設(shè)置為球域，浮心即為球域的圓心，直徑為4L。針對(duì)AUV這種不規(guī)則的幾何體，主要采用非結(jié)構(gòu)化的四面體網(wǎng)格進(jìn)行劃，為使網(wǎng)格均勻，提高阻力系數(shù)的預(yù)報(bào)精度，將靠近AUV表面流動(dòng)梯度變化較大的邊界層區(qū)域進(jìn)行加密處理。邊界層網(wǎng)格按照幾何級(jí)數(shù)膨脹的形式，采用棱柱形網(wǎng)格，其中棱柱層層數(shù)為10層，層厚設(shè)置為0.2 mm，第一層網(wǎng)格高度設(shè)置為0.05 mm,y+值為6。全域網(wǎng)格共有105萬，單元質(zhì)量平均值為0.83，相對(duì)質(zhì)量較高。網(wǎng)格劃分結(jié)果如圖6所示。

圖6 網(wǎng)格劃分結(jié)果示意圖Fig.6 Meshing result

數(shù)值計(jì)算中攻角α范圍為-20°～20°，速度v范圍 0.5～4 m/s，CFD求解器采用二階迎風(fēng)差分格式，殘差設(shè)定為10-4。通過105組仿真并將數(shù)據(jù)匯總得到阻力系數(shù)Cd隨攻角α變化曲線以及升阻比S隨攻角α變化曲線，如圖7、圖8。

由圖7可以看出阻力系數(shù)Cd隨攻角和航速的增大而增大。當(dāng)攻角α=0°時(shí)，AUV的迎流面積最小，所以此時(shí)阻力最小，迎流面積隨著攻角的增大也相應(yīng)增大，AUV的阻力同樣隨之增大。由于AUV并不是完全對(duì)稱的結(jié)構(gòu)，故阻力系數(shù)Cd的變化趨勢(shì)在α<0°和α>0°時(shí)并不是完全對(duì)稱的。

圖7 阻力系數(shù)隨攻角變化曲線Fig.7 Variation curve of drag coefficient with angle of attack

圖8 升阻比曲線Fig.8 Lift drag ratio curve

從圖8中可以看出升阻比隨攻角的增大而增大，但是升阻比最大不超過0.9，即AUV需要消耗更多的能源來維持運(yùn)動(dòng)，因此需要對(duì)其外形進(jìn)行優(yōu)化。從升阻比曲線變化趨勢(shì)中可以看出當(dāng)攻角α<-15°和α>15°時(shí)曲線增幅變緩，考慮到圖7中阻力系數(shù)的變化趨勢(shì)，選擇4 m/s、α=15°為后續(xù)優(yōu)化節(jié)點(diǎn)，此時(shí)阻力系數(shù)Cd=0.281 34，升力系數(shù)Cl=0.190 58。當(dāng)攻角α=±15°時(shí)，AUV的壓力云圖如圖9。AUV的外形有一定的對(duì)稱性，故其表面的壓力分布狀況相似，由云圖可知AUV表面出現(xiàn)了局部的壓力中心，這導(dǎo)致了誘導(dǎo)阻力的產(chǎn)生，限制了升阻比的提升。

4 優(yōu)化分析

直接采用CFD方法進(jìn)行優(yōu)化迭代計(jì)算量過大，耗時(shí)較長(zhǎng)，因此考慮構(gòu)建AUV外形結(jié)構(gòu)參數(shù)至水動(dòng)力參數(shù)的響應(yīng)面[11]，通過目標(biāo)函數(shù)計(jì)算、算法尋優(yōu)獲得最優(yōu)的水動(dòng)力外形，優(yōu)化流程如圖10所示。優(yōu)化設(shè)計(jì)的目的:調(diào)整殼體幾何尺寸，提高升阻比，改善水動(dòng)力性能，同時(shí)降低AUV的質(zhì)量。該多目標(biāo)優(yōu)化問題的數(shù)學(xué)描述如下：

目標(biāo)函數(shù)：maxSx(x)，minMx(x)

約束條件：ai≤xi≤bi

其中：設(shè)計(jì)變量X=[xi,x2,x3,x4]為殼體主要幾何參數(shù) [d3,d4,D,L]；目標(biāo)函數(shù)Sx和Mx分別為升阻比和結(jié)構(gòu)質(zhì)量；ai和bi為各幾何參數(shù)的取值上下限，即初始參數(shù)的10%。

航行器采用了魚雷狀的線型設(shè)計(jì)，故艏部段長(zhǎng)度d3和艉部段長(zhǎng)度d4以及回轉(zhuǎn)半徑D的大小直接影響了航行器艏部和艉部的形狀，中間段長(zhǎng)度L是各個(gè)功能艙長(zhǎng)度之和，對(duì)航行器的長(zhǎng)縱比起主要影響。

圖9 壓力云圖Fig.9 Pressure nephogram

圖10 優(yōu)化流程框圖Fig.10 Optimization flow chart

4.1 優(yōu)化設(shè)計(jì)方法

響應(yīng)面分析法是解決多變量問題的一種統(tǒng)計(jì)方法。為明確設(shè)計(jì)變量與響應(yīng)值之間的關(guān)系，本文中采用拉丁超立方采樣方法(LHD)[12]，采用等概率隨機(jī)正交分布的原則布置試驗(yàn)樣本點(diǎn)，以便采用較少的樣本點(diǎn)得到較高精度的響應(yīng)面模型。

響應(yīng)面擬合是在試驗(yàn)樣本點(diǎn)的基礎(chǔ)上，通過建立某些函數(shù)關(guān)系來擬合出設(shè)計(jì)變量與響應(yīng)值之間的關(guān)系[13]。選擇二階多項(xiàng)式進(jìn)行擬合，數(shù)學(xué)表達(dá)式為：

(1)

為評(píng)價(jià)響應(yīng)面的預(yù)測(cè)能力，使用決定系數(shù)R2和均方根差σ對(duì)響應(yīng)面的擬合精度進(jìn)行評(píng)估[14]，計(jì)算公式：

(2)

(3)

式中:yi為設(shè)計(jì)點(diǎn)的觀測(cè)值，當(dāng)R2趨近于1，σ趨近于0時(shí)，響應(yīng)面的擬合精度則越高。

響應(yīng)面的擬合優(yōu)度曲線如圖11所示。

圖11 擬合優(yōu)度曲線Fig.11 Goodness of fit diagram

從圖11中可以看出，不同目標(biāo)變量的響應(yīng)面預(yù)測(cè)值與其觀測(cè)值是線性變化的關(guān)系，即響應(yīng)面的擬合精度較高，擬合精度評(píng)估結(jié)果如表2所示。目標(biāo)變量的決定系數(shù)均大于0.97，均方根誤差均遠(yuǎn)小于0.1，因此響應(yīng)面的擬合精度是符合要求的。設(shè)計(jì)變量d3、d4、D、L對(duì)升力系數(shù)Cl、阻力系數(shù)Cd、質(zhì)量M和升阻比S的相對(duì)重要程度由圖12給出，中間段長(zhǎng)度L較其他3個(gè)設(shè)計(jì)變量對(duì)Cl、Cd、M和S的影響均是最大的；艏部段長(zhǎng)度d3對(duì)升、阻力系數(shù)的影響較艉部段長(zhǎng)度d4大，且艉部段長(zhǎng)度d4對(duì)升、阻力系數(shù)的作用是相反的，同時(shí)升、阻力系數(shù)受殼體回轉(zhuǎn)半徑D的影響較艉部段長(zhǎng)度d4大，且回轉(zhuǎn)半徑D對(duì)升、阻力系數(shù)的作用也是相反的；回轉(zhuǎn)半徑D對(duì)升阻比S的影響較艏部段長(zhǎng)度d3大，而艏部段長(zhǎng)度d3較艉部段長(zhǎng)度d4大；對(duì)于質(zhì)量M的影響：艏部段長(zhǎng)度d3較艉部段長(zhǎng)度d4大，而d4較回轉(zhuǎn)半徑D大。

表2 擬合精度評(píng)估結(jié)果

圖12 靈敏度直方圖Fig.12 Sensitivity diagram

圖13、圖14分別給出了升阻比S、質(zhì)量M隨艏部段長(zhǎng)度d3、艉部段長(zhǎng)度d4、殼體回轉(zhuǎn)半徑D、中間段長(zhǎng)度L的變化趨勢(shì)曲面。由圖13可以看出，升阻比隨著艏部段長(zhǎng)度d3的增加而增加，而隨著艉部段長(zhǎng)度d4的增加而降低，因此要使升阻比盡可能達(dá)到最大，d3、d4則應(yīng)在接近區(qū)間端點(diǎn)的位置取值。同時(shí)升阻比會(huì)隨著中間段長(zhǎng)度L、回轉(zhuǎn)半徑D的變化而產(chǎn)生顯著變化，L、D同樣應(yīng)該在接近區(qū)間端點(diǎn)的位置處取值，且在區(qū)間右端點(diǎn)處取值能夠?qū)崿F(xiàn)升阻比的最大化，當(dāng)4個(gè)參數(shù)共同作用時(shí)，中間段長(zhǎng)度L將對(duì)升阻比的變化起著主要作用。同理，由圖14可以看出，質(zhì)量M與4個(gè)變量呈線性變化的趨勢(shì)，且受中間段長(zhǎng)度L影響顯著，4個(gè)參變量均應(yīng)該在區(qū)間左端點(diǎn)處取值。

通過以上分析可以看出4個(gè)設(shè)計(jì)變量之間存在著互斥的關(guān)系，并考慮到實(shí)際加工制造和內(nèi)部裝配空間上的一些限制，升阻比S和質(zhì)量M無法同時(shí)達(dá)到最優(yōu)值，因此需要進(jìn)行算法尋優(yōu)，將所有可能解構(gòu)成Pareto最優(yōu)解集[16]，它是一個(gè)以犧牲其他目標(biāo)變量為代價(jià)來提高某一個(gè)目標(biāo)變量的解的集合。

圖13 輸出為S時(shí)的響應(yīng)面Fig.13 Response surface when output is S

圖14 輸出為M時(shí)的響應(yīng)面Fig.14 Response surface when output is M

4.2 多目標(biāo)遺傳算法與結(jié)果分析

為了找出目標(biāo)變量之間的最優(yōu)平衡點(diǎn)，采用偏移哈默斯利抽樣(SHS)技術(shù)選取400個(gè)樣本點(diǎn)，通過多目標(biāo)遺傳算法對(duì)優(yōu)化設(shè)計(jì)問題進(jìn)行尋優(yōu)。MOGA算法的其他參數(shù)設(shè)置如表3所示。計(jì)算升阻比S和質(zhì)量M是否滿足終止條件，若滿足則產(chǎn)生終止種群，若不滿足則通過適值分配等過程得到下一代種群，繼續(xù)判斷直到滿足終止條件[18]。通過多次迭代得到如圖15所示的Pareto最優(yōu)解集分布，圓點(diǎn)的顏色代表結(jié)果的優(yōu)劣，深藍(lán)色為最優(yōu)，黃色為最劣，其中紅色橢圓區(qū)域?yàn)樽顑?yōu)邊界，表4為擇優(yōu)選出的3個(gè)Pareto最優(yōu)解。

選取數(shù)值計(jì)算中升阻比最大，且質(zhì)量相對(duì)較小的的優(yōu)化點(diǎn)3為最終的設(shè)計(jì)結(jié)果，此時(shí)的阻力系數(shù)Cd=0.272 53，升力系數(shù)Cl=0.214 52，優(yōu)化設(shè)計(jì)前后的詳細(xì)數(shù)據(jù)如表5所示。優(yōu)化后的升阻比S和質(zhì)量M都優(yōu)于初始值，表明優(yōu)化后的AUV外形的水動(dòng)力性能優(yōu)于原外形。

圖15 Pareto最優(yōu)解集分布圖Fig.15 Distribution of Pareto solution sets

表3 MOGA參數(shù)Table 3 MOGA parameters

表4 Pareto最優(yōu)解

表5 優(yōu)化前后數(shù)據(jù)

為了驗(yàn)證優(yōu)化結(jié)果的準(zhǔn)確性，將優(yōu)化點(diǎn)中的設(shè)計(jì)變量的數(shù)據(jù)傳入FLUENT重新進(jìn)行數(shù)值計(jì)算[19]，其中，升阻比S的驗(yàn)證結(jié)果如表6所示。從表中可以看出優(yōu)化給出的結(jié)果相對(duì)于數(shù)值計(jì)算的結(jié)果最大相對(duì)誤差為1.15%，誤差較小，表明優(yōu)化結(jié)果可靠。最終制作出來的實(shí)際外形如圖16所示。

表6 升阻比的驗(yàn)證結(jié)果

圖16 實(shí)際外形圖Fig.16 Actual outline drawing

5 結(jié)論

以混合驅(qū)動(dòng)AUV為研究對(duì)象，通過CFD數(shù)值計(jì)算的方式分析了阻力系數(shù)和升阻比的變化情況，基于響應(yīng)面分析法對(duì)AUV外形水動(dòng)力性能的提升展開優(yōu)化設(shè)計(jì)，通過正交試驗(yàn)設(shè)計(jì)、響應(yīng)面擬合、靈敏度分析及MOGA算法對(duì)外形參數(shù)進(jìn)行優(yōu)化研究，結(jié)論如下：

1) 在各參數(shù)取值范圍內(nèi)，同等長(zhǎng)度條件下，取較大回轉(zhuǎn)半徑有利于提高升阻比，但也易使結(jié)構(gòu)質(zhì)量增加。AUV的長(zhǎng)度越大，阻力系數(shù)增幅明顯，不利于升阻比的提高和結(jié)構(gòu)質(zhì)量的減輕，因此應(yīng)適當(dāng)減小。在同等半徑和長(zhǎng)度條件下，增加艏部段長(zhǎng)度或減小艉部段長(zhǎng)度利于提高升阻比。

2) 基于MOGA算法搜索了全局最優(yōu)解，得到了AUV外形的優(yōu)化方案，并通過CFD方法驗(yàn)證了優(yōu)化結(jié)果的正確性。優(yōu)化后的外形較優(yōu)化前相比，升阻比由0.678提高到0.788，殼體質(zhì)量由26.3 kg降低到24.25 kg。該優(yōu)化設(shè)計(jì)方法具有良好的尋優(yōu)能力，可為AUV水動(dòng)力外形設(shè)計(jì)提供參考。