基于響應(yīng)面方法的AUV水動力外形優(yōu)化

何 隆,張 亞,李世中,王芙蓉,袁澤慧,常慧珠

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

1 引言

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

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

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

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

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

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

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

3 水動力參數(shù)計算

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

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

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

圖3 水動力模型示意圖Fig.3 Hydrodynamic model

3.1 數(shù)值計算方法驗證

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

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

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

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

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

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

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

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

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

圖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ōu)化。從升阻比曲線變化趨勢中可以看出當攻角α<-15°和α>15°時曲線增幅變緩，考慮到圖7中阻力系數(shù)的變化趨勢，選擇4 m/s、α=15°為后續(xù)優(yōu)化節(jié)點，此時阻力系數(shù)Cd=0.281 34，升力系數(shù)Cl=0.190 58。當攻角α=±15°時，AUV的壓力云圖如圖9。AUV的外形有一定的對稱性，故其表面的壓力分布狀況相似，由云圖可知AUV表面出現(xiàn)了局部的壓力中心，這導(dǎo)致了誘導(dǎo)阻力的產(chǎn)生，限制了升阻比的提升。

4 優(yōu)化分析

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

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

約束條件：ai≤xi≤bi

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

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

圖9 壓力云圖Fig.9 Pressure nephogram

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

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

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

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

(1)

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

(2)

(3)

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

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

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

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

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

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

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

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

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

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

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

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

選取數(shù)值計算中升阻比最大，且質(zhì)量相對較小的的優(yōu)化點3為最終的設(shè)計結(jié)果，此時的阻力系數(shù)Cd=0.272 53，升力系數(shù)Cl=0.214 52，優(yōu)化設(shè)計前后的詳細數(shù)據(jù)如表5所示。優(yōu)化后的升阻比S和質(zhì)量M都優(yōu)于初始值，表明優(yōu)化后的AUV外形的水動力性能優(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ōu)化結(jié)果的準確性，將優(yōu)化點中的設(shè)計變量的數(shù)據(jù)傳入FLUENT重新進行數(shù)值計算[19]，其中，升阻比S的驗證結(jié)果如表6所示。從表中可以看出優(yōu)化給出的結(jié)果相對于數(shù)值計算的結(jié)果最大相對誤差為1.15%，誤差較小，表明優(yōu)化結(jié)果可靠。最終制作出來的實際外形如圖16所示。

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

圖16 實際外形圖Fig.16 Actual outline drawing

5 結(jié)論

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

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

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