基于NFFD 和高斯近似的螺旋槳多目標(biāo)優(yōu)化

劉 旭，周喜寧，朱耀龍

(1.渤海船舶職業(yè)學(xué)院，遼寧 葫蘆島 125101；2.招商局重工(江蘇)有限公司，江蘇 南通 226000；3.華中科技大學(xué)，湖北 武漢 423201)

0 引言

螺旋槳設(shè)計(jì)及優(yōu)化一直是研究船舶快速化、經(jīng)濟(jì)性的緊迫任務(wù)之一。但由于螺旋槳參數(shù)眾多、建模繁瑣、水動(dòng)力計(jì)算周期長等問題[1]，長期限制船舶設(shè)計(jì)初期對(duì)螺旋槳的高效優(yōu)化設(shè)計(jì)。建立一套包括螺旋槳幾何變形和重構(gòu)、水動(dòng)力性能預(yù)測(cè)近似模型以及多目標(biāo)優(yōu)化的螺旋槳構(gòu)型優(yōu)化設(shè)計(jì)方法[2]，對(duì)促進(jìn)螺旋槳高效優(yōu)化設(shè)計(jì)具有重要意義。

螺旋槳的快速建模、高效優(yōu)化始終是國內(nèi)外學(xué)者關(guān)注的熱點(diǎn)研究。Ramin 等[3]采用Non-Gradinet-based Algorithm 和梯度法分別對(duì)弦長、厚度以及螺距沿半徑的分布進(jìn)行數(shù)值優(yōu)化；Florian 等[4]基于B-spline 對(duì)螺距、拱度、厚度、弦長、曲率和傾角沿半徑分布進(jìn)行擬合，并以最大效率和最小脈動(dòng)壓力為目標(biāo)進(jìn)行了螺旋槳優(yōu)化。葉禮裕[5]結(jié)合定常面元法和克里格近似理論建立螺旋槳參數(shù)與性能特征的預(yù)測(cè)模型，并利用粒子群算法對(duì)螺旋槳參數(shù)進(jìn)行了多目標(biāo)優(yōu)化設(shè)計(jì)。王超[6]用試驗(yàn)設(shè)計(jì)方法為神經(jīng)網(wǎng)絡(luò)近似模型的建立提供足量且高質(zhì)量的訓(xùn)練樣本，利用遺傳算法可以獲得更廣解集范圍更廣，而且耗時(shí)間大幅度縮減。馬丹萍[2]借助于ISIGHT 優(yōu)化平臺(tái)，集成螺旋槳的FFD 方法和CFD水動(dòng)力數(shù)值評(píng)估，提出了一套從構(gòu)型-數(shù)值求解-優(yōu)化的螺旋槳設(shè)計(jì)優(yōu)化方法。

針對(duì)目前螺旋槳研究中存在的螺旋槳模型構(gòu)造周期長、水動(dòng)力性能評(píng)估耗時(shí)多的問題，本文結(jié)合基于NURBS 基函數(shù)的自由變形算法（NURBS-based freefrom deformation，NFFD）、多輸出高斯近似模型以及NSGA-Ⅱ多目標(biāo)優(yōu)化算法，構(gòu)建一套包括螺旋槳變形重構(gòu)—水動(dòng)力性能快速預(yù)測(cè)—多目標(biāo)優(yōu)化的螺旋槳高效自動(dòng)優(yōu)化方法，即以NFFD 技術(shù)實(shí)現(xiàn)螺旋槳三維模型的變形與重構(gòu)，采用有限元數(shù)值仿真得出樣本螺旋槳的水動(dòng)力性能，基于多輸出高斯近似理論建立螺旋槳水動(dòng)力性能預(yù)測(cè)近似模型，結(jié)合近似模型和NSGA-Ⅱ算法對(duì)KP505 槳進(jìn)行多目標(biāo)優(yōu)化設(shè)計(jì)，驗(yàn)證該方法的可行性。

1 NFFD 方法及螺旋槳變形控制

基于Sederberg 等[7]提出的自由曲面變形（free form deformation，F(xiàn)FD）方法，研究人員將FFD 技術(shù)擴(kuò)展為很多不同形式[8–9]，基本包括：構(gòu)建參數(shù)控制體；將變形物嵌入?yún)?shù)控制體；參數(shù)控制體變形；計(jì)算變形物形變后的幾何坐標(biāo)[10]?；贔FD 方法，Lamousin 等[8]提出NFFD 算法，將頂點(diǎn)加權(quán)與非均勻框架結(jié)合，引入權(quán)因子和節(jié)點(diǎn)作為變量，增加了變形操作的自由度，從而達(dá)到更好的變形。馬丹萍[2]、張人會(huì)[10]、李冬琴[11]、WU[12]、馬明生[13]、唐靜[14]已將NFFD 技術(shù)應(yīng)用于船型與螺旋槳變形控制與重構(gòu)，并在船型與螺旋槳優(yōu)化設(shè)計(jì)方面取得了豐碩成果。

螺旋槳上任意一點(diǎn)絕對(duì)坐標(biāo)可以用局部坐標(biāo)系(T,U,V)表示為：

式中：Q0為局部部坐標(biāo)原點(diǎn)；(t,u,v)為點(diǎn)Q的局部坐標(biāo)；(T,U,V)為局部坐標(biāo)系標(biāo)準(zhǔn)方向向量；即控制體基本坐標(biāo)系。螺旋槳曲面參數(shù)化控制節(jié)點(diǎn)Pi,j,k可表示為：

式中：l，m，n分別表示局部坐標(biāo)在(T,U,V)3 個(gè)方向的等分份額，i=0,1,···l；j=0,1,···m；k=0,1,···n。

更進(jìn)一步，螺旋槳上任一點(diǎn)絕對(duì)坐標(biāo)可以用局部坐標(biāo)表示為：

式中：（v，u，t）為局部坐標(biāo)；Pi,j,k為參數(shù)化控制節(jié)點(diǎn)；Bil(t)，Bjm(u)，Bkn(v)為Bernstein 基函數(shù)，表示為：

NFFD 方法利用NURBS 基函數(shù)替換Bernstein 函數(shù)，螺旋槳任意一點(diǎn)Q坐標(biāo)可表示為：

式中，Wi,j,k為參數(shù)化控制節(jié)點(diǎn)Pi,j,k的權(quán)重因子。

在式（5）基礎(chǔ)上采用deBoor-Cox 遞推定義[15]和Deboor 算法[16]推導(dǎo)得出NURBS 體三次矩陣表示方法：

最后將螺旋槳嵌入以上參數(shù)控制體，并對(duì)采樣點(diǎn)進(jìn)行歸一化處理，利用牛頓法逆向搜索非線性方程組的最優(yōu)解，即可得到螺旋槳采樣點(diǎn)局部坐標(biāo)。

在建立螺旋槳與參數(shù)體的相互關(guān)系后，保持每個(gè)采樣點(diǎn)在參數(shù)體中相應(yīng)的局部坐標(biāo)(t,u,v)不變，通過改變參數(shù)化控制節(jié)點(diǎn)Pi,j,k或其權(quán)重因子Wi,j,k即可得到新的帶權(quán)控制點(diǎn)矩陣，代入式（6）中即可求解得出螺旋槳變形后的局部坐標(biāo)。

以KP505 槳為例，利用Matlab 實(shí)現(xiàn)了螺旋槳幾何形變與重構(gòu)，初始輸入螺旋槳表面離散點(diǎn)集的坐標(biāo)值，并選取0.2R，0.4R，0.6R，0.7R和0.9R處的葉寬、厚度、螺距比、曲率、傾角、最大厚度位置共30 個(gè)參數(shù)作為螺旋槳形狀控制設(shè)計(jì)變量，輸出螺旋槳變形后的表面離散點(diǎn)集的坐標(biāo)值。KP505 槳形狀控制效果和網(wǎng)格模型如圖1 所示。

圖1 KP505 槳形狀控制效果和網(wǎng)格模型Fig.1 KP505 propeller shape control effect and mesh model

2 性能預(yù)測(cè)高斯近似模型

2.1 多輸出高斯近似模型

近似模型是基于不同模型函數(shù)對(duì)設(shè)計(jì)變量和設(shè)計(jì)目標(biāo)組成的樣本數(shù)據(jù)進(jìn)行近似擬合，具有計(jì)算成本低、函數(shù)關(guān)系明確等特點(diǎn)[17]。常用的近似模型方法包括：響應(yīng)面法（RSM）、克里格法（Kringing）和人工神經(jīng)網(wǎng)絡(luò)技術(shù)[18]。高斯近似是基于Kringing 和貝葉斯定理的一種機(jī)器學(xué)習(xí)方法，對(duì)少量樣本的高階非線性問題具有良好的適應(yīng)性[19–22]。

給定訓(xùn)練樣本：

其中，任意一組樣本(xi,yi)是隨機(jī)且相互獨(dú)立的，通過近似逼近，輸出向量yi∈Y?R與輸入向量xi∈X?RP之間存在一些未知函數(shù)關(guān)系：

式中：Kn(X,X)為n×n階對(duì)稱半正定協(xié)方差矩陣，該矩陣的元素Kij=k(xp,Xq)表示xp與Xq之間的相關(guān)性；I為n×n階單位矩陣。

多輸出高斯近似模型協(xié)方差矩陣K表示為：

式中，Kf和Kx分別表示輸出和輸入向量的相關(guān)性。

訓(xùn)練樣本的輸出向量yi和測(cè)試點(diǎn)xtest對(duì)應(yīng)的目標(biāo)輸出f(xtest)服從如下聯(lián)合先驗(yàn)分布：

式中，k(xtest,X)和k(X,x)分別為n×1 和1×n 階向量。在得到y(tǒng)和f(xtest)之間的聯(lián)合先驗(yàn)概率分布之后，可以分析得到f(xtest)的后驗(yàn)概率分布如下式：

本研究中的輸入向量X和輸出向量Y分別表示為：

式中：a～f分別表示某一剖面處的葉寬、厚度、螺距比、曲率、傾角、最大厚度位置；i=0.2R，0.4R，0.6R，0.7R和0.9R，KT，10KQ，η0以及Pmin分別表示螺旋槳推力系數(shù)、扭矩系數(shù)、效率以及最小空化壓力。

2.2 敞水性能數(shù)值仿真

以STAR-CCM+13.0 流體動(dòng)力學(xué)軟件對(duì)KP505槳的敞水性能進(jìn)行仿真計(jì)算。利用NFFD 建立的原始螺旋槳曲面生成實(shí)體模型，構(gòu)建非結(jié)構(gòu)化六面體網(wǎng)格，采用k-ωSST 模型，并基于MRF 方法求解雷諾平均Navier-Stokes 方程計(jì)算不同進(jìn)速下的螺旋槳性能。

以進(jìn)速系數(shù)J=0.893 為例，在計(jì)算單元達(dá)到26 000 000時(shí)，螺旋槳的端渦和根渦現(xiàn)象如圖2 所示，即在此時(shí)實(shí)現(xiàn)了網(wǎng)格收斂。

圖2 高雷諾數(shù)下螺旋槳的端渦和根渦現(xiàn)象Fig.2 End vortex and root vortex phenomenon of propeller at high Reynolds number

圖3 為進(jìn)速系數(shù)J=0.2～1.1 時(shí)，螺旋槳推力系統(tǒng)、仿真計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比，仿真與試驗(yàn)所得推力系數(shù)KT、扭矩系數(shù)10KQ、效率η0的均方差分別為1.63%，2.41%，0.63%，說明以該軟件和網(wǎng)格數(shù)量計(jì)算的螺旋槳敞水性能精度可信。

圖3 螺旋槳仿真計(jì)算結(jié)果與試驗(yàn)結(jié)果對(duì)比Fig.3 Comparison between the simulation calculation results of the propeller and the test results

為給高斯近似模型提供初始樣本，以0.2R，0.4R，0.6R，0.7R和0.9R處的葉寬、厚度、螺距、曲率、傾角、最大厚度位置共30 個(gè)參數(shù)作為螺旋槳形狀控制設(shè)計(jì)變量，采用NFFD 技術(shù)生成100 個(gè)螺旋槳計(jì)算樣本，并計(jì)算進(jìn)速系數(shù)J=0.893 時(shí)KT，10KQ，η0以及Pmin。

2.3 近似模型的驗(yàn)證

對(duì)100 個(gè)初始樣本采用交叉驗(yàn)證方法對(duì)訓(xùn)練后的多輸出高斯代理模型進(jìn)行評(píng)估。通過交叉驗(yàn)證得到的100 個(gè)測(cè)試樣本預(yù)測(cè)誤差如圖4 所示。

圖4 高斯近似模型的交叉驗(yàn)證預(yù)測(cè)誤差Fig.4 Cross-validation prediction error of Gaussian approximation model

結(jié)果表明，各測(cè)試樣本KT，10KQ，η0以及Pmin預(yù)測(cè)值與數(shù)值模擬結(jié)果的相對(duì)誤差在3%以內(nèi)，精度滿足要求。訓(xùn)練后的多輸出高斯代理模型可用于多目標(biāo)優(yōu)化中目標(biāo)函數(shù)值的評(píng)估。

3 螺旋槳多目標(biāo)優(yōu)化

3.1 NSGA-Ⅱ多目標(biāo)優(yōu)化

NSGA-Ⅱ是由Kalyanmoy 等[23]針對(duì)NSGA 提出的一種改進(jìn)進(jìn)化算法，該算法具有多目標(biāo)適度簡化、最優(yōu)解拓展、求解精度高、高魯棒性等優(yōu)點(diǎn)，尤其針對(duì)少量樣本、多變量多目標(biāo)優(yōu)化問題就較快的求解速度。

螺旋槳敞水性能多目標(biāo)優(yōu)化問題可以描述為：

式中：葉寬0≤xa≤0.27，葉片厚度0.08 ≤xb≤0.095，螺距比1.2 ≤xc≤2.9，曲率?0.01 ≤xd≤0.02，傾角?0.5 ≤xe≤0.3，最大厚度位置0.45xb≤xf≤0.75xb，Ae螺旋槳展開面積，Ad面螺旋槳盤面積。

利用以100 個(gè)樣本建立的多輸出高斯近似模型計(jì)算螺旋槳敞水性能，并以NSGA-Ⅱ算法進(jìn)行螺旋槳的優(yōu)化設(shè)計(jì)，將種群交叉遺傳概率設(shè)置為0.92，變異概率設(shè)置為0.08。經(jīng)過1 000 代的演化，可得到Pareto 非支配最優(yōu)解集。

初始種群和Pareto 解集的扭矩系數(shù)和效率如圖5所示。可以看出最優(yōu)解是在原始設(shè)計(jì)空間之外的，其中Pareto 前沿上A點(diǎn)為綜合最優(yōu)設(shè)計(jì)，表示Pareto 前沿中的最優(yōu)綜合性能點(diǎn)，同時(shí)兼顧了扭矩系數(shù)和效率。

圖5 初始種群和Pareto 解集扭矩系數(shù)和效率Fig.5 Torque coefficient and efficiency of initial population and Pareto solution set

3.2 優(yōu)化設(shè)計(jì)結(jié)果及分析

選取優(yōu)化結(jié)果中A點(diǎn)進(jìn)行根據(jù)優(yōu)化結(jié)果分析，原模型的效率提高了2.6%、扭矩降低了3.3%。此外，盤比降低了2.5%，但螺旋槳模型的空化特性并未惡化，因?yàn)橥七M(jìn)器旋轉(zhuǎn)區(qū)域的最小壓力增加了7%。優(yōu)化方案4 槳葉r/R=0.7 處剖面輪廓如圖6 所示?？梢钥闯觯?jīng)過優(yōu)化，葉片的厚度和寬度發(fā)生了明顯的變化，在相對(duì)半徑0.7 處變小，但傾角有所增加。

圖6 優(yōu)化方案4 槳葉r/R=0.7 處剖面輪廓Fig.6 Section profile at r/R=0.7 of the optimized scheme 4 blade

槳葉寬度、剖面厚度、傾斜角沿徑向的分布如圖7所示。可以看出，經(jīng)過優(yōu)化，槳葉厚度、寬度、傾斜角在相對(duì)半徑r/R=0.7 處明顯變小，而曲率在此處所有增加。

圖7 KP505 模型槳葉寬度、厚度和傾斜角沿半徑分布Fig.7 KP505 model blade width,thickness and inclination angle distribution along the radius

優(yōu)化前后螺距比和曲率沿徑向的分布如圖8 所示。優(yōu)化前后螺旋槳的螺距比趨勢(shì)相似，而且各半徑處螺距比較對(duì)相近。此外，優(yōu)化結(jié)果表明，螺旋槳螺距比在靠近槳葉根部較大。從推進(jìn)設(shè)計(jì)的角度來看，螺旋槳螺距比在靠近槳葉根部的逐漸增大可以增加推力，從而提高效率。但是，通過增加螺距比，需要提高螺旋槳軸功率以及螺旋槳流動(dòng)區(qū)域的最小壓力。優(yōu)化前后螺旋槳曲率沿半徑分布特征相似，優(yōu)化后的螺旋槳槳葉根部曲率變?yōu)樨?fù)值，而且優(yōu)化后的曲率變小，有利于降低螺旋槳噪聲。

圖8 KP505 模型槳葉螺距比和曲率沿半徑分布Fig.8 KP505 model blade pitch ratio and curvature distribution along the radius

4 結(jié)語

1）以有限元仿真結(jié)果為樣本，以螺旋槳參數(shù)沿半徑的分布為輸入，基于多輸出高斯近似模型建立螺旋槳水動(dòng)力性能預(yù)測(cè)模型。該模型與數(shù)值仿真結(jié)果得出的螺旋槳推力系數(shù)和效率誤差在3%以內(nèi)，可以較為準(zhǔn)確預(yù)測(cè)螺旋槳水動(dòng)力性能。

2）應(yīng)用本文提出的螺旋槳高效自動(dòng)優(yōu)化方法，開展KP505 槳的優(yōu)化設(shè)計(jì)，實(shí)現(xiàn)了對(duì)KP505 槳降低扭矩系數(shù)、提高效率的多目標(biāo)優(yōu)化，形成了Pareto 解集，獲得了理論上的最佳設(shè)計(jì)方案。

3）根據(jù)對(duì)優(yōu)化結(jié)果的分析，0.65R～0.95R處螺距比和0.4R～0.6R處傾斜角對(duì)KP505 槳的效率有較大影響。

綜上所述，本文方法可以有效提高螺旋槳優(yōu)化的效率和自動(dòng)化程度，并能得到最佳設(shè)計(jì)方案。在后續(xù)工作中，將進(jìn)一步探討應(yīng)用該方法對(duì)考慮船后流場(chǎng)響應(yīng)的螺旋槳優(yōu)化設(shè)計(jì)的有效性。