孫鳴遠，劉昕瑤，張 皓，嚴 俊

(1.南方海洋科學與工程廣東省實驗室(湛江), 廣東 湛江 524000； 2.中國核動力研究設計院，成都 610213)

1 引言

螺旋槳處于船舶后的非均勻流場中，螺旋槳旋轉(zhuǎn)過程中，不同半徑處槳葉剖面的來流攻角會由于流場的不均勻性而隨時間變化，從而導致螺旋槳受到的水動力呈現(xiàn)一定的脈動特征，這是螺旋槳產(chǎn)生激振力的主要原因。螺旋槳噪聲與激振力密切相關，而在大量的試驗數(shù)據(jù)和理論研究表明，通過改善螺旋槳的側(cè)斜分布可以很大程度上降低螺旋槳產(chǎn)生的激振力。這是由于具有良好側(cè)斜分布的螺旋槳可以使不同半徑處的槳葉剖面不同時進入船后的高伴流區(qū)域，從而有效的降低螺旋槳在非均勻流場中產(chǎn)生的激振力，同時也可以降低螺旋槳由于振動等產(chǎn)生的機械噪聲[1-4]，從而減小船舶航行過程中產(chǎn)生噪聲對海洋牧場、深海養(yǎng)殖、海洋生物生態(tài)等漁業(yè)產(chǎn)業(yè)的影響。

而螺旋槳的槳葉螺距分布直接影響著各半徑槳葉剖面的迎流攻角，對螺旋槳的水動力性能有著重要影響。螺旋槳設計中要求振動和噪聲盡可能低同時螺旋槳效率較高，然而在螺旋槳實際設計過程中，降低激振力和提高螺旋槳效率往往是矛盾的，這需要在效率指標和振動指標上進行權衡優(yōu)化設計。同時隨著船舶節(jié)能減排方向的研究發(fā)展，國際拖曳水池會議(ITTC)[5-6]對于螺旋槳的優(yōu)化方向也逐漸受到關注。

螺旋槳優(yōu)化中對螺旋槳槳葉參數(shù)擬合采取的方法包括樣條曲線法[7-8]、多項式擬合[9]等方法，基本可以較為精確的對槳葉參數(shù)進行擬合。在螺旋槳數(shù)值計算方法上大部分采用面元法[10-12]或CFD方法[13-14]，仿真理論已經(jīng)較為成熟，相對面元法而言，CFD方法由于可以較為真實的模擬螺旋槳空泡形成、脫落、潰滅等現(xiàn)象，已越來越多的應用在螺旋槳優(yōu)化方向。在優(yōu)化理論的選擇上較為常用的包括模擬退火算法[15]、粒子群算法[16-17]、遺傳算法[18]等，其中NSGA II算法由于其良好的全局最優(yōu)性、種群多樣性、計算復雜度低等優(yōu)點成為多目標優(yōu)化中的經(jīng)典算法之一。

為探究七葉側(cè)斜螺旋槳非定?？张萘髦薪档图ふ窳吞岣呗菪龢实膬?yōu)化平衡問題，采用CFD數(shù)值方法結(jié)合NSGA Ⅱ多目標優(yōu)化算法對某七葉側(cè)斜螺旋槳的敞水效率和激振力進行了優(yōu)化，螺旋槳側(cè)斜及螺距沿槳葉半徑分布通過B樣條曲線進行擬合，分析結(jié)果表明所采用優(yōu)化策略對提高七葉側(cè)斜螺旋槳敞水效率和降低激振力是有效的。

2 螺旋槳優(yōu)化方法

2.1 基于精英策略的非支配排序遺傳算法

遺傳算法在求解非確定性多項式問題上相比傳統(tǒng)的搜索算法有明顯的優(yōu)勢。非支配排序遺傳算法[19](non-dominated sorting genetic algorithm，NSGA)通過樣本中個體之間的支配關系分層實現(xiàn)，求解得到的非劣最優(yōu)解的分布較為均勻，但是計算過程比較復雜且無精英策略。改進的基于精英策略的非支配排序遺傳算法NSGA Ⅱ[20](non-dominated sorting genetic algorithm Ⅱ，NSGA Ⅱ)對NSGA進行了改進，基于精英策略的算法大幅改進了程序計算復雜、共享參數(shù)依賴性較大等問題，大幅提高了計算效率。同時NSGA Ⅱ算法在選擇交叉的個體時，省去了共享函數(shù)、共享半徑等的求解過程，實現(xiàn)快速非支配排序。同時NSGA Ⅱ算法引入了擁擠距離的概念用來保證最優(yōu)解集的多樣性和均勻性，結(jié)合其精英策略保證種群進化過程中不會丟失最優(yōu)解，NSGA Ⅱ算法相比于NSGA算法有很大的提升。

2.2 螺旋槳非定常水動力性能預報方法

螺旋槳非定常水動力性能預報采用CFD數(shù)值模擬方法，對三維不可壓縮流動，控制方程如下：

(1)

對流場中流體微元的動量守恒方程可以寫為：

(2)

式中:ρ為流體密度；u、v、w為流體微元的速度分量；p為流體微元壓力；Fx、Fy、Fz分別為流體微元三個方向的質(zhì)量力。

采用SSTk-ω兩方程湍流模型，該模型結(jié)合了標準k-ω方程和標準k-ε方程的優(yōu)點，在模擬螺旋槳旋轉(zhuǎn)等復雜流動時可以在近場和遠場模擬中都取得較好的效果，相關研究在文獻[21-22]中有詳細論述，在此不詳細說明。

2.3 螺旋槳幾何參數(shù)化表達及優(yōu)化理論

對螺旋槳進行多目標優(yōu)化首先需要對螺旋槳的幾何模型進行參數(shù)化表達，常用的參數(shù)化表達方法包括B樣條曲線擬合[7-8,23-24]、貝塞爾函數(shù)曲線擬合[25]、多項式擬合[26]等。B樣條曲線擬合方法由于有著局部可控性、曲線光順性等優(yōu)勢從而有著廣泛的應用。本文采用B樣條函數(shù)對螺旋槳側(cè)斜及螺距沿槳葉半徑的分布進行曲線擬合，并將B樣條曲線的控制點作為設計變量。圖1所示為七葉側(cè)斜螺旋槳曲線擬合分布與原始分布。圖中螺距分布由N1～N4四個點進行控制，側(cè)斜分布由N5～N9五個點進行控制。從圖中可以看出，所選取B樣條曲線控制點個數(shù)擬合螺旋槳幾何參數(shù)分布的精度是可以接受的。

圖.1 螺旋槳參數(shù)擬合與原參數(shù)分布擬合曲線Fig.1 Comparison of propeller parameter fitting and original parameter distribution

NSGA優(yōu)化算法最初由Srinivas提出，其理論核心包括基于非支配排序原理的種群分級和基于共享小生境原理的虛擬擁擠度計算2個方面，在多目標優(yōu)化領域有著廣泛應用，但隨著復雜度的增加，NAGA算法存在計算復雜度過高、需要人為制定共享參數(shù)等缺陷。因此，在前人的研究基礎上，Deb等學者提出了改進型的NAGA-Ⅱ型算法，提出了精英保留、非支配排序等概念，同時提出了擁擠度的概念，保證了種群的多樣性，大大提高了算法的適用度，同時也擴展了Pareto解集的分布范圍。

擁擠度是指種群中給定個體的周圍個體密度，直觀上可表示為個體i周圍僅僅包含個體i本身時的最大長方形的長，如圖2所示，優(yōu)化目標函數(shù)為f1和f2，個體i的擁擠度為d。擁擠度的引入實際上是為了保證在交叉運算中不至于出現(xiàn)基因過度集中的現(xiàn)象。

圖2 擁擠度概念示意圖Fig.2 Schematic diagram of congestion

兩目標優(yōu)化問題的擁擠距離的計算方法如下：

(3)

式中：L[i]d為第i個個體的第m個目標函數(shù)值；L為個體之間的距離。

2.4 螺旋槳側(cè)斜及螺距優(yōu)化方法

所進行的螺旋槳側(cè)斜及螺距優(yōu)化，是基于2.2節(jié)中螺旋槳非定常水動力數(shù)值計算方法及2.1節(jié)中多目標優(yōu)化算法，將2.3節(jié)中的擬合參數(shù)(N1～N9)在合理范圍內(nèi)取值，從而得到一系列滿足設計目標及約束條件的螺旋槳設計方案。該優(yōu)化方法的具體做法如下：

優(yōu)化變量：

N1～N9

優(yōu)化目標：

min(KFX(N1,N2,…,N9))

(4)

min(-η0(N1,N2,…,N9))

(5)

約束條件：

(6)

因此螺旋槳優(yōu)化的流程為：通過CFD數(shù)值方法預報得到螺旋槳敞水推進效率以及軸向一倍頻非定常力系數(shù)等優(yōu)化目標，由NSGA Ⅱ算法進行優(yōu)化計算，根據(jù)優(yōu)化目標及約束條件滿足情況，重新設計優(yōu)化變量，得到新的螺旋槳設計方案再次開始CFD數(shù)值計算，得到優(yōu)化目標后再次進行優(yōu)化結(jié)果評判、重新設計優(yōu)化變量，直至滿足優(yōu)化迭代要求為止。

3 數(shù)值方法及計算模型

3.1 數(shù)值方法驗證

螺旋槳優(yōu)化選取原型槳為某七葉側(cè)斜螺旋槳E1619，螺旋槳幾何參數(shù)如表1所示，幾何模型如圖3所示。

表1 E1619槳幾何參數(shù)

圖3 E1619槳幾何模型示意圖Fig.3 E1619 propeller geometric model

螺旋槳非定常水動力數(shù)值計算域為圓柱體，計算中將計算域劃分為靜止域和旋轉(zhuǎn)域2部分，旋轉(zhuǎn)域同樣為圓柱體，計算域尺寸及邊界條件劃分如圖4所示。靜止域與旋轉(zhuǎn)域之間邊界條件為Interface，旋轉(zhuǎn)域與靜止域之間的數(shù)值傳遞通過Interface進行，因此為保證數(shù)據(jù)傳遞的準確性，需保證Interface的網(wǎng)格大小基本一致，同時在Interface兩側(cè)各設置一層網(wǎng)格大小相同的邊界層，網(wǎng)格劃分結(jié)果如圖5所示。

圖4 E1619槳計算域劃分示意圖Fig.4 E1619 propeller computational domain division

圖5 網(wǎng)格劃分結(jié)果示意圖Fig.5 Mesh division results

圖6 數(shù)值計算結(jié)果與試驗結(jié)果曲線Fig.6 Comparison of numerical calculation results and experimental results

3.2 網(wǎng)格收斂性分析

網(wǎng)格劃分情況對計算結(jié)果的影響較大，因此數(shù)值計算需要對網(wǎng)格進行收斂性分析，只有當網(wǎng)格數(shù)量的增加對計算結(jié)果影響很小時，此時的計算結(jié)果才具有可信度。參照文獻[26-27]的方法進行網(wǎng)格收斂性分析，不同網(wǎng)格方案滿足一定細化率要求：

k=(Nf/Nc)1/d

(7)

其中，k為細化率，取為1.2；Nf為加密后網(wǎng)格數(shù)量；Nc為加密前網(wǎng)格數(shù)量；d為網(wǎng)格空間維度，取為3。通過上述網(wǎng)格細化策略得到3組不同的網(wǎng)格方案Mesh a、Mesh b和Mesh c，三組網(wǎng)格數(shù)量分別為6.53×106,1.16×107,2.10×107。網(wǎng)格劃分基本尺寸均為螺旋槳直徑D，槳葉表面網(wǎng)格邊界層均設置為6層，網(wǎng)格劃分參數(shù)如表2所示。

表2 網(wǎng)格劃分參數(shù)

按照上述網(wǎng)格劃分方案，計算螺旋槳設計工況(J=0.7)時的水動力特性，表3給出了3種網(wǎng)格方案計算得到的螺旋槳水動力參數(shù)。可以看出，隨著網(wǎng)格的加密，數(shù)值計算結(jié)果與試驗值的差異逐漸減小。采用Mesh b網(wǎng)格時，推力系數(shù)與試驗值誤差為2.16%，敞水效率與試驗值誤差為0.7%，已經(jīng)滿足計算精度要求，因此后續(xù)計算中均采用Mesh b網(wǎng)格進行，槳葉表面網(wǎng)格劃分如圖7所示。

表3 不同網(wǎng)格劃分下KT and η0與試驗值參數(shù)

圖7 螺旋槳表面網(wǎng)格劃分示意圖Fig.7 Mesh division of the propeller surface

4 計算結(jié)果及討論

對N1～N9共計9個優(yōu)化變量在螺旋槳設計工況(J=0.7)處對螺旋槳進行多目標優(yōu)化，優(yōu)化目標如式(4)、式(5)，約束條件ε=5%，即螺旋槳推力系數(shù)KT范圍為KT≥0.232 8。9個優(yōu)化變量在±10%范圍內(nèi)進行變化，設定種群大小為15，迭代次數(shù)為20次，共生成300個個體，其中滿足設定約束條件的可行解有232個。圖8為可行解的氣泡狀分布圖，圖中氣泡直徑表示螺旋槳扭矩系數(shù)，氣泡顏色表示螺旋槳推力系數(shù)。圖中紅色虛線圈中顏色為紅色及粉紅色圓圈表示不滿足推力系數(shù)約束條件的解，剩余部分為Pareto解，可以看出Pareto解集在Pareto前沿上分布較為均勻。

圖8 可行解氣泡狀分布圖Fig.8 Bubble diagram distribution of feasible solutions

參照圖8計算結(jié)果在Pareto前沿上選取四型螺旋槳進行分析：P236、P177、P102和P42。其中P236為滿足推力系數(shù)約束條件下一階軸向脈動力指標最優(yōu)的方案；P177為滿足推力系數(shù)約束條件下敞水效率性能最優(yōu)的方案；P102為滿足推力系數(shù)約束條件下綜合考慮一階軸向脈動力指標及敞水效率性能的折中方案；P42為不考慮推力系數(shù)約束條件下的最優(yōu)方案。表4為上述幾型螺旋槳設計參數(shù)及設計目標參數(shù)。

表4 優(yōu)化槳與原槳幾何參數(shù)

從表4中可以看出：優(yōu)化螺旋槳敞水效率均較原槳有一定程度提高，一階軸向力脈動值均較原槳下降，同時都相對原槳有一定程度的推力損失。在考慮推力系數(shù)約束條件限制下，P236槳敞水效率提高0.8%，一階軸向力脈動水平下降10%；P177槳敞水效率提高1.7%，一階軸向力脈動水平下降2.8%；折中方案P102槳敞水效率提高1.4%，一階軸向力脈動水平下降6.7%。P42槳敞水效率提高2.4%，一階軸向力脈動水平下降9.6%，是所有優(yōu)化方案中優(yōu)化效果最好的，可惜不滿足推力系數(shù)限制條件，因此該方案不可取。

圖9和圖10分別為3組優(yōu)化槳和原槳的螺距比分布和側(cè)斜分布曲線。

圖9 優(yōu)化槳與原槳螺距分布曲線Fig.9 Comparison of the pitch distribution between the optimized propeller and the original propeller

從圖9中可以看出，除個別半徑處優(yōu)化槳螺距比原槳略大外，優(yōu)化槳整體螺距比分布相較原槳總體呈減小的趨勢。其中P236槳和P177槳0.7R處螺距相比原槳降低了9.9%和11.3%，折中方案P102槳較原槳降低了4.8%。

圖10 優(yōu)化槳與原槳側(cè)斜分布曲線Fig.10 Comparison of the skew distribution between the optimized propeller and the original propeller

從圖10可以看出，在r大于0.7R后滿足推力系數(shù)約束條件的3組優(yōu)化槳的側(cè)斜均較原槳有所提高。P236槳和P177槳0.9R處的槳葉側(cè)斜角度較原槳增加4.2%和4.8%，折中方案P102槳較原槳增加了2.2%。整體來看優(yōu)化槳與原槳的側(cè)斜分布變化幅度較小。上述數(shù)據(jù)在一定程度上說明，將原型槳降低各半徑處螺距比、提高槳葉靠近葉梢位置處的側(cè)斜角度是提高推進效率、降低一階軸向脈動力水平有益的探索方向。

圖11為原槳與3組優(yōu)化槳推力系數(shù)時歷曲線，推力系數(shù)在一個穩(wěn)定值附近隨時間脈動，推力系數(shù)時歷曲線符合螺旋槳推力脈動特性。從圖中已標出推力系數(shù)脈動幅值來看，3組優(yōu)化槳中P236槳脈動幅值最小，P177槳脈動幅值最大，P102槳則介于上述二者之間，該結(jié)果與3組優(yōu)化槳一階軸向力脈動計算結(jié)果相符。

圖11 原槳和優(yōu)化槳的推力系數(shù)時歷曲線Fig.11 Comparison of thrust coefficient time-history curves of the original propeller and the optimized propeller

5 結(jié)論

1) 采用基于CFD數(shù)值方法與多目標優(yōu)化方法結(jié)合的方式，對某七葉側(cè)斜螺旋槳進行螺距及側(cè)斜分布等多目標全局優(yōu)化， Pareto解集在Pareto前沿上分布均勻，表明所采用的NAGS-Ⅱ多目標優(yōu)化方法可有效提高螺旋槳推進效率，降低一階軸向力，所采取的優(yōu)化策略可行。

2) 優(yōu)化后的螺旋槳總體螺距相比原型槳較小，側(cè)斜分布在r大于0.7R后有所增加，在考慮螺旋槳激振力更優(yōu)情況下選擇P236槳，一階軸向力脈動值下降10%，考慮效率更優(yōu)時選擇P177槳，敞水效率可提高1.7%。

3) 七葉側(cè)斜螺旋槳優(yōu)化中，將原型槳降低各半徑處螺距比(r<0.7R)、提高槳葉靠近葉梢位置處(r>0.7R)的側(cè)斜角度是提高推進效率、降低一階軸向脈動力的研究方向。