基于NFFD 和高斯近似的螺旋槳多目標優化

劉 旭，周喜寧，朱耀龍

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

0 引言

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

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

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

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

基于Sederberg 等[7]提出的自由曲面變形（free form deformation，FFD）方法，研究人員將FFD 技術擴展為很多不同形式[8–9]，基本包括：構建參數控制體；將變形物嵌入參數控制體；參數控制體變形；計算變形物形變后的幾何坐標[10]。基于FFD 方法，Lamousin 等[8]提出NFFD 算法，將頂點加權與非均勻框架結合，引入權因子和節點作為變量，增加了變形操作的自由度，從而達到更好的變形。馬丹萍[2]、張人會[10]、李冬琴[11]、WU[12]、馬明生[13]、唐靜[14]已將NFFD 技術應用于船型與螺旋槳變形控制與重構，并在船型與螺旋槳優化設計方面取得了豐碩成果。

螺旋槳上任意一點絕對坐標可以用局部坐標系(T,U,V)表示為：

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

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

更進一步，螺旋槳上任一點絕對坐標可以用局部坐標表示為：

式中：（v，u，t）為局部坐標；Pi,j,k為參數化控制節點；Bil(t)，Bjm(u)，Bkn(v)為Bernstein 基函數，表示為：

NFFD 方法利用NURBS 基函數替換Bernstein 函數，螺旋槳任意一點Q坐標可表示為：

式中，Wi,j,k為參數化控制節點Pi,j,k的權重因子。

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

最后將螺旋槳嵌入以上參數控制體，并對采樣點進行歸一化處理，利用牛頓法逆向搜索非線性方程組的最優解，即可得到螺旋槳采樣點局部坐標。

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

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

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

2 性能預測高斯近似模型

2.1 多輸出高斯近似模型

近似模型是基于不同模型函數對設計變量和設計目標組成的樣本數據進行近似擬合，具有計算成本低、函數關系明確等特點[17]。常用的近似模型方法包括：響應面法（RSM）、克里格法（Kringing）和人工神經網絡技術[18]。高斯近似是基于Kringing 和貝葉斯定理的一種機器學習方法，對少量樣本的高階非線性問題具有良好的適應性[19–22]。

給定訓練樣本：

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

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

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

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

訓練樣本的輸出向量yi和測試點xtest對應的目標輸出f(xtest)服從如下聯合先驗分布：

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

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

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

2.2 敞水性能數值仿真

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

以進速系數J=0.893 為例，在計算單元達到26 000 000時，螺旋槳的端渦和根渦現象如圖2 所示，即在此時實現了網格收斂。

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

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

圖3 螺旋槳仿真計算結果與試驗結果對比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 個參數作為螺旋槳形狀控制設計變量，采用NFFD 技術生成100 個螺旋槳計算樣本，并計算進速系數J=0.893 時KT，10KQ，η0以及Pmin。

2.3 近似模型的驗證

對100 個初始樣本采用交叉驗證方法對訓練后的多輸出高斯代理模型進行評估。通過交叉驗證得到的100 個測試樣本預測誤差如圖4 所示。

圖4 高斯近似模型的交叉驗證預測誤差Fig.4 Cross-validation prediction error of Gaussian approximation model

結果表明，各測試樣本KT，10KQ，η0以及Pmin預測值與數值模擬結果的相對誤差在3%以內，精度滿足要求。訓練后的多輸出高斯代理模型可用于多目標優化中目標函數值的評估。

3 螺旋槳多目標優化

3.1 NSGA-Ⅱ多目標優化

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

螺旋槳敞水性能多目標優化問題可以描述為：

式中：葉寬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 個樣本建立的多輸出高斯近似模型計算螺旋槳敞水性能，并以NSGA-Ⅱ算法進行螺旋槳的優化設計，將種群交叉遺傳概率設置為0.92，變異概率設置為0.08。經過1 000 代的演化，可得到Pareto 非支配最優解集。

初始種群和Pareto 解集的扭矩系數和效率如圖5所示。可以看出最優解是在原始設計空間之外的，其中Pareto 前沿上A點為綜合最優設計，表示Pareto 前沿中的最優綜合性能點，同時兼顧了扭矩系數和效率。

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

3.2 優化設計結果及分析

選取優化結果中A點進行根據優化結果分析，原模型的效率提高了2.6%、扭矩降低了3.3%。此外，盤比降低了2.5%，但螺旋槳模型的空化特性并未惡化，因為推進器旋轉區域的最小壓力增加了7%。優化方案4 槳葉r/R=0.7 處剖面輪廓如圖6 所示。可以看出，經過優化，葉片的厚度和寬度發生了明顯的變化，在相對半徑0.7 處變小，但傾角有所增加。

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

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

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

優化前后螺距比和曲率沿徑向的分布如圖8 所示。優化前后螺旋槳的螺距比趨勢相似，而且各半徑處螺距比較對相近。此外，優化結果表明，螺旋槳螺距比在靠近槳葉根部較大。從推進設計的角度來看，螺旋槳螺距比在靠近槳葉根部的逐漸增大可以增加推力，從而提高效率。但是，通過增加螺距比，需要提高螺旋槳軸功率以及螺旋槳流動區域的最小壓力。優化前后螺旋槳曲率沿半徑分布特征相似，優化后的螺旋槳槳葉根部曲率變為負值，而且優化后的曲率變小，有利于降低螺旋槳噪聲。

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

4 結語

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

2）應用本文提出的螺旋槳高效自動優化方法，開展KP505 槳的優化設計，實現了對KP505 槳降低扭矩系數、提高效率的多目標優化，形成了Pareto 解集，獲得了理論上的最佳設計方案。

3）根據對優化結果的分析，0.65R～0.95R處螺距比和0.4R～0.6R處傾斜角對KP505 槳的效率有較大影響。

綜上所述，本文方法可以有效提高螺旋槳優化的效率和自動化程度，并能得到最佳設計方案。在后續工作中，將進一步探討應用該方法對考慮船后流場響應的螺旋槳優化設計的有效性。