孫鳴遠,劉昕瑤,張 皓,嚴 俊
(1.南方海洋科學與工程廣東省實驗室(湛江), 廣東 湛江 524000; 2.中國核動力研究設計院,成都 610213)
螺旋槳處于船舶后的非均勻流場中,螺旋槳旋轉過程中,不同半徑處槳葉剖面的來流攻角會由于流場的不均勻性而隨時間變化,從而導致螺旋槳受到的水動力呈現一定的脈動特征,這是螺旋槳產生激振力的主要原因。螺旋槳噪聲與激振力密切相關,而在大量的試驗數據和理論研究表明,通過改善螺旋槳的側斜分布可以很大程度上降低螺旋槳產生的激振力。這是由于具有良好側斜分布的螺旋槳可以使不同半徑處的槳葉剖面不同時進入船后的高伴流區域,從而有效的降低螺旋槳在非均勻流場中產生的激振力,同時也可以降低螺旋槳由于振動等產生的機械噪聲[1-4],從而減小船舶航行過程中產生噪聲對海洋牧場、深海養殖、海洋生物生態等漁業產業的影響。
而螺旋槳的槳葉螺距分布直接影響著各半徑槳葉剖面的迎流攻角,對螺旋槳的水動力性能有著重要影響。螺旋槳設計中要求振動和噪聲盡可能低同時螺旋槳效率較高,然而在螺旋槳實際設計過程中,降低激振力和提高螺旋槳效率往往是矛盾的,這需要在效率指標和振動指標上進行權衡優化設計。同時隨著船舶節能減排方向的研究發展,國際拖曳水池會議(ITTC)[5-6]對于螺旋槳的優化方向也逐漸受到關注。
螺旋槳優化中對螺旋槳槳葉參數擬合采取的方法包括樣條曲線法[7-8]、多項式擬合[9]等方法,基本可以較為精確的對槳葉參數進行擬合。在螺旋槳數值計算方法上大部分采用面元法[10-12]或CFD方法[13-14],仿真理論已經較為成熟,相對面元法而言,CFD方法由于可以較為真實的模擬螺旋槳空泡形成、脫落、潰滅等現象,已越來越多的應用在螺旋槳優化方向。在優化理論的選擇上較為常用的包括模擬退火算法[15]、粒子群算法[16-17]、遺傳算法[18]等,其中NSGA II算法由于其良好的全局最優性、種群多樣性、計算復雜度低等優點成為多目標優化中的經典算法之一。
為探究七葉側斜螺旋槳非定常空泡流中降低激振力和提高螺旋槳效率的優化平衡問題,采用CFD數值方法結合NSGA Ⅱ多目標優化算法對某七葉側斜螺旋槳的敞水效率和激振力進行了優化,螺旋槳側斜及螺距沿槳葉半徑分布通過B樣條曲線進行擬合,分析結果表明所采用優化策略對提高七葉側斜螺旋槳敞水效率和降低激振力是有效的。
遺傳算法在求解非確定性多項式問題上相比傳統的搜索算法有明顯的優勢。非支配排序遺傳算法[19](non-dominated sorting genetic algorithm,NSGA)通過樣本中個體之間的支配關系分層實現,求解得到的非劣最優解的分布較為均勻,但是計算過程比較復雜且無精英策略。改進的基于精英策略的非支配排序遺傳算法NSGA Ⅱ[20](non-dominated sorting genetic algorithm Ⅱ,NSGA Ⅱ)對NSGA進行了改進,基于精英策略的算法大幅改進了程序計算復雜、共享參數依賴性較大等問題,大幅提高了計算效率。同時NSGA Ⅱ算法在選擇交叉的個體時,省去了共享函數、共享半徑等的求解過程,實現快速非支配排序。同時NSGA Ⅱ算法引入了擁擠距離的概念用來保證最優解集的多樣性和均勻性,結合其精英策略保證種群進化過程中不會丟失最優解,NSGA Ⅱ算法相比于NSGA算法有很大的提升。
螺旋槳非定常水動力性能預報采用CFD數值模擬方法,對三維不可壓縮流動,控制方程如下:
(1)
對流場中流體微元的動量守恒方程可以寫為:
(2)
式中:ρ為流體密度;u、v、w為流體微元的速度分量;p為流體微元壓力;Fx、Fy、Fz分別為流體微元三個方向的質量力。
采用SSTk-ω兩方程湍流模型,該模型結合了標準k-ω方程和標準k-ε方程的優點,在模擬螺旋槳旋轉等復雜流動時可以在近場和遠場模擬中都取得較好的效果,相關研究在文獻[21-22]中有詳細論述,在此不詳細說明。
對螺旋槳進行多目標優化首先需要對螺旋槳的幾何模型進行參數化表達,常用的參數化表達方法包括B樣條曲線擬合[7-8,23-24]、貝塞爾函數曲線擬合[25]、多項式擬合[26]等。B樣條曲線擬合方法由于有著局部可控性、曲線光順性等優勢從而有著廣泛的應用。本文采用B樣條函數對螺旋槳側斜及螺距沿槳葉半徑的分布進行曲線擬合,并將B樣條曲線的控制點作為設計變量。圖1所示為七葉側斜螺旋槳曲線擬合分布與原始分布。圖中螺距分布由N1~N4四個點進行控制,側斜分布由N5~N9五個點進行控制。從圖中可以看出,所選取B樣條曲線控制點個數擬合螺旋槳幾何參數分布的精度是可以接受的。

圖.1 螺旋槳參數擬合與原參數分布擬合曲線Fig.1 Comparison of propeller parameter fitting and original parameter distribution
NSGA優化算法最初由Srinivas提出,其理論核心包括基于非支配排序原理的種群分級和基于共享小生境原理的虛擬擁擠度計算2個方面,在多目標優化領域有著廣泛應用,但隨著復雜度的增加,NAGA算法存在計算復雜度過高、需要人為制定共享參數等缺陷。因此,在前人的研究基礎上,Deb等學者提出了改進型的NAGA-Ⅱ型算法,提出了精英保留、非支配排序等概念,同時提出了擁擠度的概念,保證了種群的多樣性,大大提高了算法的適用度,同時也擴展了Pareto解集的分布范圍。
擁擠度是指種群中給定個體的周圍個體密度,直觀上可表示為個體i周圍僅僅包含個體i本身時的最大長方形的長,如圖2所示,優化目標函數為f1和f2,個體i的擁擠度為d。擁擠度的引入實際上是為了保證在交叉運算中不至于出現基因過度集中的現象。

圖2 擁擠度概念示意圖Fig.2 Schematic diagram of congestion
兩目標優化問題的擁擠距離的計算方法如下:

(3)
式中:L[i]d為第i個個體的第m個目標函數值;L為個體之間的距離。
所進行的螺旋槳側斜及螺距優化,是基于2.2節中螺旋槳非定常水動力數值計算方法及2.1節中多目標優化算法,將2.3節中的擬合參數(N1~N9)在合理范圍內取值,從而得到一系列滿足設計目標及約束條件的螺旋槳設計方案。該優化方法的具體做法如下:
優化變量:
N1~N9
優化目標:
min(KFX(N1,N2,…,N9))
(4)
min(-η0(N1,N2,…,N9))
(5)
約束條件:
(6)

因此螺旋槳優化的流程為:通過CFD數值方法預報得到螺旋槳敞水推進效率以及軸向一倍頻非定常力系數等優化目標,由NSGA Ⅱ算法進行優化計算,根據優化目標及約束條件滿足情況,重新設計優化變量,得到新的螺旋槳設計方案再次開始CFD數值計算,得到優化目標后再次進行優化結果評判、重新設計優化變量,直至滿足優化迭代要求為止。
螺旋槳優化選取原型槳為某七葉側斜螺旋槳E1619,螺旋槳幾何參數如表1所示,幾何模型如圖3所示。

表1 E1619槳幾何參數

圖3 E1619槳幾何模型示意圖Fig.3 E1619 propeller geometric model
螺旋槳非定常水動力數值計算域為圓柱體,計算中將計算域劃分為靜止域和旋轉域2部分,旋轉域同樣為圓柱體,計算域尺寸及邊界條件劃分如圖4所示。靜止域與旋轉域之間邊界條件為Interface,旋轉域與靜止域之間的數值傳遞通過Interface進行,因此為保證數據傳遞的準確性,需保證Interface的網格大小基本一致,同時在Interface兩側各設置一層網格大小相同的邊界層,網格劃分結果如圖5所示。

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

圖5 網格劃分結果示意圖Fig.5 Mesh division results


圖6 數值計算結果與試驗結果曲線Fig.6 Comparison of numerical calculation results and experimental results
網格劃分情況對計算結果的影響較大,因此數值計算需要對網格進行收斂性分析,只有當網格數量的增加對計算結果影響很小時,此時的計算結果才具有可信度。參照文獻[26-27]的方法進行網格收斂性分析,不同網格方案滿足一定細化率要求:
k=(Nf/Nc)1/d
(7)
其中,k為細化率,取為1.2;Nf為加密后網格數量;Nc為加密前網格數量;d為網格空間維度,取為3。通過上述網格細化策略得到3組不同的網格方案Mesh a、Mesh b和Mesh c,三組網格數量分別為6.53×106,1.16×107,2.10×107。網格劃分基本尺寸均為螺旋槳直徑D,槳葉表面網格邊界層均設置為6層,網格劃分參數如表2所示。

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

表3 不同網格劃分下KT and η0與試驗值參數

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

圖8 可行解氣泡狀分布圖Fig.8 Bubble diagram distribution of feasible solutions
參照圖8計算結果在Pareto前沿上選取四型螺旋槳進行分析:P236、P177、P102和P42。其中P236為滿足推力系數約束條件下一階軸向脈動力指標最優的方案;P177為滿足推力系數約束條件下敞水效率性能最優的方案;P102為滿足推力系數約束條件下綜合考慮一階軸向脈動力指標及敞水效率性能的折中方案;P42為不考慮推力系數約束條件下的最優方案。表4為上述幾型螺旋槳設計參數及設計目標參數。

表4 優化槳與原槳幾何參數
從表4中可以看出:優化螺旋槳敞水效率均較原槳有一定程度提高,一階軸向力脈動值均較原槳下降,同時都相對原槳有一定程度的推力損失。在考慮推力系數約束條件限制下,P236槳敞水效率提高0.8%,一階軸向力脈動水平下降10%;P177槳敞水效率提高1.7%,一階軸向力脈動水平下降2.8%;折中方案P102槳敞水效率提高1.4%,一階軸向力脈動水平下降6.7%。P42槳敞水效率提高2.4%,一階軸向力脈動水平下降9.6%,是所有優化方案中優化效果最好的,可惜不滿足推力系數限制條件,因此該方案不可取。
圖9和圖10分別為3組優化槳和原槳的螺距比分布和側斜分布曲線。

圖9 優化槳與原槳螺距分布曲線Fig.9 Comparison of the pitch distribution between the optimized propeller and the original propeller
從圖9中可以看出,除個別半徑處優化槳螺距比原槳略大外,優化槳整體螺距比分布相較原槳總體呈減小的趨勢。其中P236槳和P177槳0.7R處螺距相比原槳降低了9.9%和11.3%,折中方案P102槳較原槳降低了4.8%。

圖10 優化槳與原槳側斜分布曲線Fig.10 Comparison of the skew distribution between the optimized propeller and the original propeller
從圖10可以看出,在r大于0.7R后滿足推力系數約束條件的3組優化槳的側斜均較原槳有所提高。P236槳和P177槳0.9R處的槳葉側斜角度較原槳增加4.2%和4.8%,折中方案P102槳較原槳增加了2.2%。整體來看優化槳與原槳的側斜分布變化幅度較小。上述數據在一定程度上說明,將原型槳降低各半徑處螺距比、提高槳葉靠近葉梢位置處的側斜角度是提高推進效率、降低一階軸向脈動力水平有益的探索方向。
圖11為原槳與3組優化槳推力系數時歷曲線,推力系數在一個穩定值附近隨時間脈動,推力系數時歷曲線符合螺旋槳推力脈動特性。從圖中已標出推力系數脈動幅值來看,3組優化槳中P236槳脈動幅值最小,P177槳脈動幅值最大,P102槳則介于上述二者之間,該結果與3組優化槳一階軸向力脈動計算結果相符。

圖11 原槳和優化槳的推力系數時歷曲線Fig.11 Comparison of thrust coefficient time-history curves of the original propeller and the optimized propeller
1) 采用基于CFD數值方法與多目標優化方法結合的方式,對某七葉側斜螺旋槳進行螺距及側斜分布等多目標全局優化, Pareto解集在Pareto前沿上分布均勻,表明所采用的NAGS-Ⅱ多目標優化方法可有效提高螺旋槳推進效率,降低一階軸向力,所采取的優化策略可行。
2) 優化后的螺旋槳總體螺距相比原型槳較小,側斜分布在r大于0.7R后有所增加,在考慮螺旋槳激振力更優情況下選擇P236槳,一階軸向力脈動值下降10%,考慮效率更優時選擇P177槳,敞水效率可提高1.7%。
3) 七葉側斜螺旋槳優化中,將原型槳降低各半徑處螺距比(r<0.7R)、提高槳葉靠近葉梢位置處(r>0.7R)的側斜角度是提高推進效率、降低一階軸向脈動力的研究方向。