權(quán)曉波，王 惠，魏海鵬，程少華

（北京宇航系統(tǒng)工程研究所，北京100076）

基于代理模型的水下航行體頭型優(yōu)化設(shè)計方法研究

權(quán)曉波，王 惠，魏海鵬，程少華

（北京宇航系統(tǒng)工程研究所，北京100076）

文章采用Kriging代理模型技術(shù)針對水下航行體頭型優(yōu)化設(shè)計開展研究。采用CFD方法分析了不同頭型水下航行體的流體動力性能，將性能參數(shù)作為初始樣本建立代理模型；根據(jù)測試樣本點的預(yù)測均方誤差選擇加點策略更新代理模型，提高代理模型預(yù)測精度；采用非劣分類遺傳算法對阻力和表面壓力開展多目標優(yōu)化，尋找近似最優(yōu)解。結(jié)果表明：基于代理模型技術(shù)的水下航行體頭型多目標優(yōu)化設(shè)計方法可以有效提高設(shè)計效率，獲得具有良好水動特性的航行體頭型。

水下航行體；頭型設(shè)計；Kriging代理模型；優(yōu)化設(shè)計

0 引 言

頭型設(shè)計是水下航行體外形設(shè)計的主要內(nèi)容，直接決定航行體的阻力和抗空化能力，對航行體的水下運動特性產(chǎn)生重要影響[1]。阻力減小可以提高航速，增大航程，增加攜帶有效載荷的能力[2]；抗空化能力的提高可以降低流場空化形成的可能性，避免由空化產(chǎn)生的沖擊、振動及噪聲。為使水下航行體獲得更小的阻力和更高的抗空化能力，需采用多目標優(yōu)化方法對水下航行體頭型進行尋優(yōu)。

在頭型優(yōu)化設(shè)計過程中，需反復(fù)修改外形，進行多次流體動力特性計算[3]。傳統(tǒng)的頭型設(shè)計一般采用數(shù)值CFD計算方式，可獲得相對準確的計算結(jié)果，但CFD分析消耗大量的計算資源和時間，導(dǎo)致迭代優(yōu)化設(shè)計時效率較低。應(yīng)用代理模型技術(shù)可有效地解決這一問題，通過采用近似方法實現(xiàn)計算效率和計算精度的有效匹配，在優(yōu)化設(shè)計領(lǐng)域獲得廣泛的應(yīng)用。宋保維等[4]基于代理模型進行了魚雷外形阻力的單目標優(yōu)化，證明此優(yōu)化方法可以克服計算成本高的缺點，獲得較為滿意的魚雷外形；姚拴寶等[5]基于改進的蟻群算法和改進的Kriging代理模型開展了列車頭型的三維氣動減阻優(yōu)化設(shè)計研究，減小了頭部阻力。

本文針對水下航行體頭型多目標優(yōu)化設(shè)計問題，應(yīng)用代理模型技術(shù)作為學(xué)科分析工具和優(yōu)化算法之間的接口，獲取水下航行體優(yōu)化頭型。首先，基于Kriging近似方法和最優(yōu)拉丁超立方試驗設(shè)計方法建立代理模型；然后，根據(jù)最小化均方誤差方法選擇加點策略，更新代理模型以提高預(yù)測精度；而后，選取NSGA-Ⅱ遺傳算法進行頭型多目標優(yōu)化設(shè)計，獲取兼顧阻力及抗空化能力的Pareto前沿解。最終，對優(yōu)化后頭型應(yīng)用CFD數(shù)值仿真計算方法進行仿真，驗證了本方法的有效性。

1 研究方法

1.1 Kriging代理模型

代理模型技術(shù)是指采用近似模型建立輸入到輸出的響應(yīng)系統(tǒng)，包括近似方法和試驗設(shè)計方法的選取。

1.1.1 Kriging近似方法

Kriging代理模型是一種估計方差最小的無偏估計模型[6]。對于輸入x，其響應(yīng)y可表示為：

式中：β是回歸系數(shù)，f（ x）T是x的n階多項式，隨機分布的誤差為z，具有如下統(tǒng)計特性：

式中：xi，xj為樣本中的任意兩點；R（ θ, xi,xj）是以θ為參數(shù)的相關(guān)模型，表示樣本的空間相關(guān)性。本文選取的高斯相關(guān)函數(shù)表達式為：

式中：xik，xjk是樣本點xi，xj第k維元素。

對于m維初始樣本點，其響應(yīng)為Y，可得多項式F和相關(guān)矩陣R。對于一個待測點x的，可得多項表達式f和相關(guān)陣r，其響應(yīng)值利用樣本點xi的響應(yīng)值Y線性加權(quán)疊加計算，即=C（ x）TY。在無偏性要求下，模型預(yù)測期望方差表示為：

拉格朗日方程求解上述方程取最小值時可得：

C（ x）與待測點和已知點間的相關(guān)性矩陣有關(guān)。由最小二乘估計可得：

方差的極大似然估計為σ2：

1.1.2 試驗設(shè)計方法

試驗設(shè)計是構(gòu)建代理模型的基礎(chǔ)，通過試驗設(shè)計給出一組具有代表性的小量樣本設(shè)計點，盡可能地反映設(shè)計空間的特性[7]。例如，對于二因子6水平的研究，全因子設(shè)計（FFD）需要6×6個點，而最優(yōu)拉丁超立方設(shè)計（Opt LHD）只需研究6個點，如圖1所示?？梢钥闯觯疚倪x取的最優(yōu)拉丁超立方設(shè)計可使試驗點均勻地分布在設(shè)計空間，空間填充性和均衡性高，擬合非線性響應(yīng)效果好。

1.1.3 基于最小化均方誤差的Kriging改進模型

為了在小樣本的條件下提高Kriging近似方法的擬合精度，本文基于最小化均方誤差（MSE）的方法評估代理模型的預(yù)測效果，由原樣本點生成新增樣本點，增加樣本信息，擬合這兩部分樣本點共同組成的樣本空間，達到提高了擬合精度的目標。

首先要針對不同的輸入?yún)?shù)，進行全因子設(shè)計，假設(shè)共有m維輸入設(shè)計參數(shù)，每個輸入設(shè)計參數(shù)選取的水平數(shù)是ni，則可以獲得測試樣本的容量為n1×n2…×nm；然后應(yīng)用Kriging方法對測試樣本進行預(yù)測，分析測試樣本MSE的變化規(guī)律，找到不滿足預(yù)測要求的一組“差”點，作為補充樣本點更新原樣本點；最終，以新樣本點為基礎(chǔ)建立代理模型，驗證精度是否滿足要求，流程如圖2所示。

圖1 全因子設(shè)計和最優(yōu)拉丁超立方設(shè)計試驗點分布Fig.1 Sample distribution of FFD and Opt LHD

圖2 改進的Kriging模型的構(gòu)建流程Fig.2 Flow chart of improved Kriging model

1.2 數(shù)值分析方法

本文采用基于結(jié)構(gòu)網(wǎng)格和雷諾平均N-S方程求解航行體阻力系數(shù)和表面壓力系數(shù)。具體計算方法簡要描述為：使用FLUENT軟件求解三維定常不可壓縮流動，湍流模型為k-ω的SST形式，壓力速度耦合方式采用SIMPLEC算法，對流項采用2階迎風(fēng)格式?？紤]流場對稱性，采用ICEM軟件對一半流場劃分結(jié)構(gòu)化網(wǎng)格。在靠近旋成體表面對網(wǎng)格加密，形成邊界層，最靠近物面一層網(wǎng)格間距為10-4m量級。

為了驗證本文數(shù)值仿真方法的正確性，對某水下航行體表面壓力系數(shù)進行數(shù)值模擬，網(wǎng)格和壓力云圖如圖3所示。航行體表面壓力系數(shù)的仿真結(jié)果與試驗結(jié)果對比如圖4所示，由圖中可以看出仿真結(jié)果與試驗結(jié)果相差較小，壓力系數(shù)極小值的誤差僅為2.9%，表明本文采用的CFD計算方法具有較高的準確度。

從壓力云圖可以看出，頭部迎背水面存在低壓區(qū)，易產(chǎn)生空化現(xiàn)象。以空化初生時空化數(shù)的大小描述航行體產(chǎn)生空化現(xiàn)象的難易程度：初生空化數(shù)越小，說明該物體的抗空化能力越高。不同物體在液體中運動，初生空化數(shù)由最小壓力系數(shù)cpmin的絕對值確定[8]，越小，初生空化數(shù)越小，該物體抗空化能力越高。因此，本文中選取航行體表面最小壓力系數(shù)來表征水下航行體頭型抗空化能力的大小，通過尋找最小壓力系數(shù)絕對值的最小值，獲取抗空化能力最強的航行體頭型。

圖3 水下航行體結(jié)構(gòu)化網(wǎng)格和壓力云圖Fig.3 Structured mesh and pressure contour of underwater vehicle

圖4 航行體表面壓力系數(shù)仿真結(jié)果與試驗結(jié)果對比Fig.4 Comparison of prediction value and CFD value of Cp

1.3 優(yōu)化算法

遺傳算法（GA）是一種借鑒生物自然選擇和遺傳機制的隨機搜索方法[9]，由于具有群體搜索特性并且不需要輔助的梯度信息，使得其具有較好的全局搜索性能，在當(dāng)今工程計算領(lǐng)域的多目標優(yōu)化問題中有著廣泛的應(yīng)用。NSGA-Ⅱ的優(yōu)點在于具有獲得最優(yōu)Pareto前沿解的能力[10]，而且Pareto前沿前進能力增強，探索性能良好，因此，本文選取NSGA-Ⅱ作為多目標優(yōu)化算法。

2 算例建模

2.1 研究對象

以某水下航行體為研究對象，外形如圖5所示。為了方便、準確地描述水下航行體的頭型曲線，需要進行頭型曲線的參數(shù)化設(shè)計，以魚雷線形擬合研究為基礎(chǔ)，選用常見的平端頭曲線的雙參數(shù)立方多項式形式作為水下航行體頭部曲線擬合的基本解析表達式：

式中：R為水下航行體半徑；X0為頭型曲線過渡段長度；k1和k2為曲線的控制參數(shù)，k1是前端點的曲率變化率，當(dāng)k1取值大時，頭部線型豐滿，反之消瘦；k2是后端點的曲率變化率，當(dāng)k2取值大時，頭部線型消瘦，反之豐滿，規(guī)律如圖6所示。

圖5 水下航行體外形Fig.5 Configuration of underwater vehicle

圖6 曲線控制參數(shù)k1、k2對頭型的影響Fig.6 Influence of curve parameters

2.2 試驗設(shè)計

設(shè)計變量選取頭型曲線過渡段長度X0∈[0.3 m,1 m]，頭型曲線的控制參數(shù)k1∈[0,1]、k2∈[0,15]，采用最優(yōu)拉丁超立方試驗設(shè)計方法獲得一組初始樣本點，樣本點容量為12個，空間分布情況如圖7所示。

2.3 仿真計算

在5°攻角、來流速度15 m/s條件下，數(shù)值計算12組初始樣本點的阻力系數(shù)與壓力系數(shù)，結(jié)果如表1。

2.4 基于Kriging改進方法建立代理模型

根據(jù)表1中12組初始樣本點，采用Kriging方法建立初始代理模型。為了提高初始代理模型的預(yù)測精度，采用最小化均方誤差的方法更新初始Kriging模型，分為以下幾步：

（1）采用全因子試驗設(shè)計方法，將設(shè)計變量k1、k2、X0分別間隔0.1、1、0.05取值，獲得11×16×15維訓(xùn)練樣本；

（2）應(yīng)用初始代理模型預(yù)測訓(xùn)練樣本的流體力系數(shù)，提取表面壓力系數(shù)cp的預(yù)測MSE，以MSE值最大的8組樣本點為補充樣本點。針對補充樣本點開展仿真，結(jié)果如表2所示；

（3）以表1的初始點和表2的補充點作為樣本點建立代理模型，獲得表面壓力系數(shù)cp的預(yù)測MSE如圖8所示?？梢钥闯?，補充8組角點數(shù)據(jù)后，角點處MSE值明顯降低，代理模型精度水平有所提升；

（4）在樣本點之外隨機生成2組測試樣本點，采用代理模型開展預(yù)測，評估代理模型的預(yù)測精度。流體力系數(shù)的預(yù)測值與仿真值誤差對比如表3所示，可以看出，阻力系數(shù)的預(yù)測誤差不超過2%，表面壓力系數(shù)誤差不超過13%，最大預(yù)測方差滿足收斂要求，可利用此代理模型代替真實物理模型進行優(yōu)化設(shè)計。

圖7 初始樣本點空間分布Fig.7 Space distribution of sample points

表1 初始樣本點阻力系數(shù)與壓力系數(shù)Tab.1 Drag coefficient and maximum pressure coefficient of design variables

表2 8組補充樣本點的阻力系數(shù)與壓力系數(shù)Tab.2 Drag and pressure coefficients of added design variables

圖8 測試樣本點的MSE值Fig.8 MSE of test point

表3 代理模型預(yù)測誤差與初始誤差對比Tab.3 Prediction error companion of improved model and origin model

3 多目標優(yōu)化結(jié)果

在來流速度v=15 m/s、攻角α=5°的條件下，水下航行體頭型優(yōu)化的多目標為：mincd、min。

采用NSGA-Ⅱ算法對上述問題進行多目標優(yōu)化，圖9給出了優(yōu)化結(jié)果的Pareto前沿圖。Pareto前沿的兩個端點分別為壓力系數(shù)絕對值的最小值和阻力系數(shù)cd的最小值，中間點則是以犧牲一個優(yōu)化目標為代價改善另一個優(yōu)化目標。將Pareto優(yōu)化前沿點與初始及補充樣本點相比可以看出，其阻力系數(shù)和表面壓力系數(shù)絕對值均較小，達到了阻力更低和抗空化能力更高的多目標優(yōu)化目的。但由于兩個優(yōu)化目標在一定范圍內(nèi)相互制約，均獲得最優(yōu)解的理想狀態(tài)不可能達到。因此，Pareto優(yōu)化前沿上的點均可作為可行優(yōu)化解。

從Pareto優(yōu)化解中選取一組特征頭型曲線，分析曲線的共同特征可知，優(yōu)化后k1較小，k2、X0值較大，頭型具有如圖10所示特征。

圖9 多目標優(yōu)化結(jié)果Fig.9 Multi-optimization result(Pareto)

圖10 優(yōu)化后頭型曲線集Fig.10 A set of optimal forehead configurations

在Pareto優(yōu)化解集中，針對壓力系數(shù)絕對值最小點開展仿真，驗證基于代理模型優(yōu)化設(shè)計的準確性，相應(yīng)的仿真計算值和優(yōu)化預(yù)測值對比如表4所示?？梢钥闯觯枇ο禂?shù)和壓力系數(shù)的仿真值與預(yù)測值差距較小，誤差分別為0.8%、-2.6%，表明本文中采用的基于代理模型的多目標優(yōu)化設(shè)計方法準確度較高，可以將該方法獲得的優(yōu)化解集作為水下航行體頭型優(yōu)化結(jié)果。

表4 Pareto邊界點頭型參數(shù)的數(shù)值仿真和優(yōu)化預(yù)測對比Tab.4 Prediction value and CFD value comparison of Pareto boundary point

圖11 阻力系數(shù)和壓力系數(shù)隨迭代次數(shù)變化規(guī)律Fig.11 Drag coefficient and maximum pressure coefficient changing with iterative steps

本文進行優(yōu)化設(shè)計時，初始人口數(shù)設(shè)置為12，遺傳代數(shù)為100，迭代次數(shù)共1 200次，阻力系數(shù)和壓力系數(shù)隨迭代次數(shù)的變化如圖11所示?？梢钥闯?，在此過程中，共調(diào)用代理模型1 200次，如果要直接使用CFD方法進行優(yōu)化設(shè)計，需要進行上千次計算，而使用代理模型代替數(shù)值仿真，調(diào)用代理模型的時間與CFD仿真時間相比幾乎可以不計。由此可見，引入代理模型可以減少優(yōu)化過程時間，實現(xiàn)水下航行體頭型優(yōu)化設(shè)計效率的大幅提升。

4 結(jié) 論

本文將改進的Kriging代理模型技術(shù)應(yīng)用于水下航行體頭型優(yōu)化設(shè)計，得到了如下結(jié)論：

（1）基于最優(yōu)拉丁超立方試驗設(shè)計及最小化均方誤差加點策略的改進Kriging代理模型可以使用數(shù)量較少的樣本點建立可靠的近似模型。與常規(guī)Kriging代理模型技術(shù)相比提高了預(yù)測MSE水平，滿足了優(yōu)化過程的收斂要求；

（2）多目標優(yōu)化可以給出同時改善阻力和抗空化能力兩方面目標的優(yōu)化前沿解集，得到了具有良好水動特性的航行體外形，與單目標優(yōu)化相比具有更好的適用性；

（3）基于代理模型技術(shù)的水下航行體頭型優(yōu)化設(shè)計可以大大減少CFD仿真計算次數(shù)，降低成本，提高設(shè)計效率，同時，優(yōu)化結(jié)果的預(yù)測值與仿真值相差較小，證明此方法是一種高效可靠的水下航行體頭型設(shè)計方法。

[1]張宇文.魚雷外形設(shè)計（上冊）[M].西安:西北工業(yè)大學(xué)出版社,1998,8:317-372. Zhang Yuwen.Torpedo shape design[M].Xi’an:Northwestern Polytechnic University Press,1998,8:317-372.

[2]王 鵬,宋保維,劉旭琳,余德海,姜 軍.魚雷外形多學(xué)科設(shè)計優(yōu)化與仿真[J].系統(tǒng)仿真學(xué)報,2008,20(7):1915-1918. Wang Peng,Song Baowei,Liu Xulin,Yu Dehai,Jiang Jun.Multidisciplinary design optimization and simulation validation for shape design of torpedo[J].Journal of System Simulation,2008,20(7):1915-1918.

[3]馬 洋,王丹丹,楊 濤,豐志偉,張青斌.基于代理模型的運載器頭罩外形優(yōu)化設(shè)計[J].航空動力學(xué)報,2014,29(1): 192-198. Ma Yang,Wang Dandan,Yang Tao,Feng Zhiwei,Zhang Qinbin.Optimization design of shroud configuration of transporter based on surrogate model[J].Journal of Aerospace Power,2014,29(1):192-198.

[4]溫慶國,宋保維,王鵬,王司令.改進的Kriging近似方法及其在魚雷外形優(yōu)化中的應(yīng)用[J].上海交通大學(xué)學(xué)報,2013, 47(4):613-618. Wen Qingguo,Song Baowei,Wang Peng,Wang Siling.An improved Kriging approximation and its application to shape design optimization of torpedo[J].Journal of Shanghai Jiao Tong University,2013,47(4):613-618.

[5]姚拴寶,郭迪龍,楊國偉.基于徑向基函數(shù)網(wǎng)格變形的高速列車頭型優(yōu)化[J].力學(xué)學(xué)報,2013,45(6):982-986. Yao Shuanbao,Guo Dilong,Yang Guowei.Aerodynamic optimization of high-speed train based on RBF mesh deformation [J].Chinese Journal of Theoretical and Applied Mechanics,2013,45(6):982-986.

[6]Shawn E.Gano,Harold Kim,Don E.Brown II.Comparison of three surrogate modeling techniques:datascape,Kriging, and Second order regression[C].11th AIAA/ISSMO Multidisciplinary Analysis Conference,6-8 Sep,2006.AIAA-2006-7048.

[7]Grosso A.Finding maximum Latin hypercube designs by iterated local search heuristics[J].J Oper Res,2009,197(2): 541-547.

[8]王獻孚.空化泡和超空化泡流動理論及應(yīng)用[M].北京:國防工業(yè)出版社,2009:3.

[9]Pratab A,Deb K.A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization:NSGA-II[R]. Kan-GAL report 200001,Indian Institute of Technology,Kanpur,India,2000.

[10]賴宇陽.Isight參數(shù)優(yōu)化理論與實例詳解[M].北京:北京航空航天大學(xué)出版社,2012:264. Lai Yuyang.Isight parameter optimization theory and example explanation[M].Beijing:Beijing University of Aeronautics and Astronautics Press,2012:264.

Configuration optimization design of underwater vehicle based on surrogate model

QUAN Xiao-bo,WANG Hui,WEI Hai-peng,CHENG Shao-hua
(Beijing Institute of Astronautical System Engineering,Beijing 100076,China)

Based on Kriging surrogate model,a configuration optimization design method of underwater vehicle was studied.A set of different forehead configurations were chosen as original sample points to construct the surrogate model.The hydrodynamic forces are calculated by CFD analysis tool.To improve the attainable precision,sampling strategy was used according to the MSE(mean square error)of test points.After that,NSGA-Ⅱoptimization method was implemented under given condition to find the approximate optimal solution considering minimum drag coefficient and maximum pressure coefficient.The results show that the multi-object optimization design method based on the surrogate model gives a set of Pareto optimal solutions.Therefore it can improve the configuration design efficiency of underwater vehicle,and find the vehicle configuration with good hydrodynamic performance.

underwater vehicle;forehead configuration;Kriging surrogate model;optimal design

TV131.2

：A

10.3969/j.issn.1007-7294.2016.10.006

1007-7294（2016）10-1262-07

2016-07-18

權(quán)曉波（1976-），男，博士后，研究員；王 惠（1993-），女，碩士研究生，E-mail:wanghui0912wh@126.com。