漁船球鼻艏參數優化研究

李廣年，于 群，杜 林，秦海瑞，袁文鑫

（1.寧波大學海運學院，浙江寧波 315211；2.上海市東方海事工程技術有限公司，上海 200011）

0 引 言

近年來，由于過度捕撈及環境污染導致近海漁業資源匱乏。為緩解我國漁業壓力、建設新型現代化漁業體系、推進海洋強國戰略部署，發展遠洋漁業勢在必行。2020年“中國遠洋漁業履約白皮書”提出支持漁船更新改造，鼓勵采用安全高效、節能減排的技術和船用裝備。船舶參數化水動力構型為綠色節能船型的研發提供了極大的便利，是船舶行業實現“碳達峰”、“碳中和”目標的技術途徑之一。球鼻艏可在一定程度上減小船舶阻力，但其泛化性較弱，即同一球鼻艏應用到其他船時并不一定能夠獲得滿意的減阻效果。因此，球鼻艏形狀、參數等與船體有恰當的配合是球鼻艏方案設計、優化的重點方向。

國內外研究人員針對球鼻艏做過一系列的研究。Chrismianto 等[1]使用4 個設計參數對球鼻艏進行參數化設計，利用多目標遺傳算法對KCS（KRISO Container Ship）進行優化。Lu等[2]通過參數化方法對球鼻艏進行變形，并利用改進非支配排序遺傳算法（Non-dominated Sorting Genetic Algorithms，NSGA-II）對一艘大型集裝箱船進行總阻力優化，結果表明優化后的球鼻艏使船體總阻力降低2.845%。Huang 等[3]利用徑向基函數（Radial Basis Function，RBF）和多目標人工蜂群優化算法對一艘貨船球鼻艏進行優化，并通過計算和試驗驗證得到了最優船型。Hong等[4]結合CFD（Computational Fluid Dynamics）技術和DoE（Design of Experiment）技術，以最小阻力為目標對某漁船球鼻艏進行了優化設計。Yu等[5]使用粒子群算法對一艘散貨船船艏進行優化研究，結果表明最優船型的興波阻力降低13.2%。隨著多目標與多學科優化的興起以及設計變量維度的增加，優化設計往往需要海量的計算，基于近似模型的優化技術是解決這種復雜優化的方法之一。Lin 等[6]以Kriging 模型作為替代模型，使用NSGA-II對某雙艉漁船阻力性能進行優化。Guerrero 等[7]以兩個球鼻艏參數為設計變量，使用Kriging 模型作為替代模型并進行優化，最終使阻力降低約7%。Feng 等[8]采用支持向量回歸（Support Vector Regression,SVR）替代模型和NSGA-II方法對某近海養殖船進行優化，優化船型的阻力比原始船型減小了1.59%。Coppedè 等[9]基于幾何變形方法、CFD 計算、普通Kriging 模型以及高斯回歸的響應面法搭建的優化框架對KCS進行優化，使總阻力減小4%。綜上所述，盡管專家學者們針對球鼻艏優化做了大量的工作，但探究多個球鼻艏特征參數對漁船阻力的影響并將參數化方法運用到漁船球鼻艏優化的并不多見。

本文以一艘遠洋漁船為研究對象，以設計航速（Fr=0.28）時總阻力最小為目標，介紹了一種漁船球鼻艏參數化優化方法。首先，選取漁船球鼻艏特征參數和特征曲線對其進行參數化建模，以球鼻艏的幾何參數作為設計變量，利用敏感性分析方法結合CFD技術對各變量進行靈敏度分析，甄別出對阻力性能變化最為敏感的若干變量；然后，針對敏感變量使用窮舉搜索法對樣本空間進行采樣，并通過數值仿真計算樣本總阻力；最后，利用多項式構建替代模型得到設計變量與總阻力系數之間的函數關系并使用遺傳算法進行優化，以獲得阻力最優的球鼻艏參數。

1 基本理論與方法

本文將數值預報、靈敏度分析、最優化技術等方法應用于漁船球鼻艏優化問題。其中，數值預報是建立漁船球鼻艏優化問題數學模型的基礎，數值計算精度直接影響了優化結果的質量；靈敏度分析可以去除球鼻艏模型中對阻力變化影響較小的參數，降低整個優化的計算成本及時間成本；最優化技術是解決漁船球鼻艏優化問題的必要手段，通過優化算法能夠在設計空間內快速、準確地搜尋最優解，即最佳漁船球鼻艏。

采用Numeca公司的FINE/Marine 模塊進行數值模擬，在等溫條件下，針對所求解的不可壓縮黏性流，質量和動量守恒方程分別為

式中，t為時間，ρ為密度，V為控制體，S為圍成控制體的面積，Ud為S上n方向的速度，U和p分別為速度和壓力，Ui為在xi坐標軸方向上的平均速度分量，τij和gi分別為黏性應力張量和重力矢量，Ii和Ij為方向向量，ci為i流體的體積份額。

使用k-ω(SST-Menter)湍流模型，其湍流能k和湍流耗散頻率ω滿足：

式中：xj為坐標軸；Uj為在xj坐標軸方向上的平均速度分量；μ為分子黏度；μt為湍流渦黏度；τtij為湍流雷諾應力張量；Sij為平均應變率張量；F1為輔助混合函數；Pω為ω的導出項；β*、σk和σω2分別為湍流模型常數。

在設計初期，可以通過靈敏度分析快速有效地確定靈敏度較大的參數，為后續球鼻艏優化工作節約時間及計算成本。一種典型的靈敏度分析方法是OAT(one-at-a-time)，它可以直觀地觀察輸入因素對輸出的影響。本文以球鼻艏的特征參數為設計變量，使用OAT 進行靈敏度分析，具體內容為：改變一個設計變量，其它設計變量保持不變，記錄設計變量的變化對總阻力系數的影響；然后將該設計變量返回初始值，對其它設計變量進行重復操作；最后通過分析每個設計變量對總阻力系數Ct的影響，將靈敏度最高的設計變量作為優化的設計變量。

選擇有限個點分析設計變量與目標方程的響應，即確定替代模型?？山铺娲O計變量與目標方程的關系，從而減少計算量。建立替代模型的前提是對樣本空間的選擇，窮舉搜索法（Exhaustive Search，ES）是一種具有規則結構，能夠系統地列舉多狀態網格的所有可能情況，并檢查每個情況是否符合問題描述的方法。本文利用ES方法將設計變量的約束區間細分為給定數量的等距間隔，逐一匹配確定樣本點，并通過CFD計算對應的Ct。

采用二元多項式構建替代模型，其多項式擬合模型定義如下：

式中：f(x,y)即為總阻力系數Ct；x、y為自變量；a1、a2、…、ai、…、aj為擬合系數，a0為截距項；m、n為擬合次數。替代模型生成之后，需要對其進行精度驗證，確保其可以精確地反應出設計變量對Ct的響應。擬合優度（goodness of fit）是指回歸曲線對真實值的擬合程度，度量擬合優度的統計量是確定系數Rsquare。

式中，yi為樣本的真實值，y^i為替代模型預測值，yˉi為樣本均值，SSR為預測值與真實值均值差的平方和，SSE為真實值與預測值差的平方和，SST為真實值與真實值均值差的平方和。Rsquare給出了目標變量的可變性程度，它介于0 到1 之間。Rsquare越高，說明模型擬合得越好。在擬合模型過程中，當模型Rsquare≥0.95時，目標變量對總阻力系數的響應滿足精度要求。

遺傳算法（Genetic Algorithm，GA）是借鑒生物進化過程而提出的一種啟發式搜索算法。其主要特點是直接對結構對象進行操作，且不存在求導和對函數連續性的限定，采用概率化的尋優方法，不需要確定的規則就能自動獲取和指導優化的搜索空間，自適應地調整搜索方向。其通用表達式為

式中：f(xi)為目標函數；x1,x2,x3…xn為自變量；s.t.為約束條件。本文在使用GA進行優化時，首先，以總阻力最小為目標，隨機生成目標的初代解作為初代種群；其次，對種群中的個體，即設計變量進行編碼，以替代模型作為個體的適應度，計算種群中每個個體的適應度；接著，根據適應度的高低選擇參與繁衍的父體與母體，并對被選出的父體與母體執行遺傳操作，即復制父體與母體的基因，并采用交叉、變異等算子產生出子代；最后，找出所有子代中適應度最高個體，即總阻力最小的設計變量作為最優解。

2 問題描述

在進行優化前，首先要進行CFD 精度計算和參數化建模。選擇KCS 船型進行CFD 精度驗證，數值模擬計算域尺寸及邊界條件如圖1所示。采用全六面體非結構網格，網格的生成方式是由體到面，并將物體表面附近網格適當細化并投影到物面上形成貼體網格。為了準確模擬流動和船體阻力，在自由液面處和壁面處進行網格加密，整個計算域的網格數量約為134萬，船體網格劃分如圖2所示。

圖1 邊界條件Fig.1 Simulation domain and boundary conditions

圖2 網格劃分Fig.2 Mesh division

通過對比KCS在不同傅汝德數下總阻力系數計算值與MOERI拖曳水池試驗數據，如圖3所示，可以看出二者相對吻合，平均誤差在1.13%左右。說明本文所采用的數值分析方法適用于船舶阻力計算，其計算精度可達到預期效果。

圖3 KCS阻力數據對比Fig.3 Resistance curve of KCS

本文研究對象為一艘金槍魚船，如圖4 所示，該船球鼻艏最大橫剖面位于首柱附近，球鼻艏的長度Lb為2.43 m，球鼻艏的最大橫剖面處的最大寬度Bb為1.07 m，球鼻艏縱剖面最前端點的高度Zb為1.88 m，球鼻艏中縱剖面面積Abl為5.63 m2，球鼻艏中最大橫剖面面積Abt為5.14 m2。

圖4 球鼻艏示意圖Fig.4 Profile of bulbous bow

金槍魚船的縮尺比λ=16，主要船型參數如表1所示。在球鼻艏建模過程中，針對球鼻艏主要參數及幾何特征，共有7個球鼻艏特征參數，其初始值及取值范圍如表2所示。表2中，X1和X2是與Abl相關的特征參數，X4與Abt相關。同時還設置了特征曲線，它們分別為最大寬度曲線以及最大寬度高度曲線，即Bb在船長方向上，從最大橫剖面處至球鼻艏最遠點寬度和高度的變化曲線。

表1 金槍魚船的主要船型參數Tab.1 Principal hull form parameters of the tuna seiner

表2 球鼻艏特征參數Tab.2 The original seven parameters of bulbous bow

本文使用OAT方法，在不考慮多個變量耦合的情況下對每個變量進行分析，在每個參數所在的取值區間內均分取樣，每個參數都得到5 個樣本，其樣本取值隨著編號的增大而增大。如圖5 所示，在變化區間內隨著每個設計變量的改變，其對應Ct的變化趨勢和變化范圍都有所不同。每個設計變量對應的Ct曲線起點不同，但最終近似匯集于兩處。第一處相關變量為X6和X7，X7對應的Ct不減反增，X6對應的Ct變化是波動的；X6和X7對應的Ct總體變化量分別為0.957%和-1.451%。第二處相關變量為X1～X5，各自對應的Ct總體變化量分別為-3.650%、-2.317%、-4.035%、-1.363%和-1.663%；X3對應的Ct在變化區間內單調下降，X4、X5對應的Ct變化趨勢大致相同，都是先增加再減小，X1和X2對應的Ct變化也呈下降趨勢，但在第三個樣本處略有增加。

圖5 球鼻艏特征參數靈敏度分析Fig.5 Sensitivity analysis of seven parameters

綜上所述，在不考慮變量之間的耦合情況下，對每個設計變量進行靈敏度分析。就整體而言，X1和X3對應的Ct總體變化量比其它五個參數要大得多。進一步考慮X1和X3的耦合效果對總阻力性能的影響，其他變量均取初始值。以X1和X3為優化問題的設計變量，使用ES 取樣，共生成100 個樣本，經過ES抽樣的空間點和平面投影點分布如圖6所示。

圖6 ES取樣分布Fig.6 Uniform sampling on the design space

可以看出，在X1和X3耦合的情況下，其樣本點的變化趨勢與單一變量靈敏度分析大致相同，Ct隨X1和X3增大而減小，尤其是在0.65～0.80 范圍內，Ct下降的趨勢較大。綜合考慮模型精度與計算量的情況下，最終選擇二元五次多項式作為本文的代理模型，即

其中擬合系數如表3所示。

表3 擬合系數Tab.3 Fitting coefficients

在保證精度滿足要求情況下，剔除了樣本集中個別數據奇異點，經過計算，最終替代模型的擬合優度Rsquare為0.95，擬合曲面如圖7所示。

圖7 擬合曲面Fig.7 Polynomial regression surface

從圖7 可以看出，隨著X1和X3增大，Ct呈減小的趨勢較為明顯，在f( 0.6,0.6 )附近出現極大值，最小值則在f( 0.9,0.9 )附近。

引入替代模型后，優化問題可表示為

式中，目標函數為在設計航速（Fr=0.28）時Ct最小，約束條件為|Δo-Δb|≤1%Δb（其中，Δb為母型船排水量，Δo為優化船排水量）及對應變量的取值范圍。

設置種群數量為50，迭代次數為100，交叉率和突變率分別為0.8 和0.01。GA 的進化終止條件為達到最大進化代數或進化過程中連續10代種群之間的適應度平均值變化小于設定值1E-6，優化結果如圖8所示。

圖8 GA迭代過程Fig.8 Optimization convergent curve

從圖8（a）可以看出，優化在40代以后開始收斂，且種群適應度平均值的變化小于1E-6，從圖8（b）可以看出，第40代之前個體的平均距離呈振蕩狀態，第40代以后個體的平均距離逐漸收斂，個體間的平均距離為0。

3 結果與討論

運用GA 對替代模型進行優化得到的最優解見表4，最優解所對應的球鼻艏形狀如圖9 所示。從球鼻艏的橫剖面可以看出，在橫向上其形狀變化不大，球鼻艏的最大寬度不變；在縱向上球鼻艏上半部分的形狀較母型更加飽滿，但是球鼻艏的長度沒有變化，最遠端的高度卻有所下降。

表4 優化結果Tab.4 Optimization results

圖9 球鼻艏形狀對比Fig.9 Comparison of profiles

根據表5 阻力成分（Fr= 0.28）對比可知，在球鼻艏優化過程中，與Re相關的粘性阻力沒有下降，反而增加0.16%，但與Fr相關的壓阻力下降了5.43%。

表5 阻力成分（Fr=0.28）Tab.5 Occupation changes of total resistance（Fr=0.28）

圖10 為母型船與優化船在Fr為0.20～0.30 工況下的總阻力系數對比。圖中顯示，在Fr為0.22、0.24、0.30 時總阻力下降較多，分別減小了3.32%、3.27%和3.57%；而在Fr為0.20、0.26、0.28時，總阻力下降相對較小，分別下降1.35%、2.96%和2.68%。可以看出，本文優化后的球鼻艏在低航速時（Fr≤0.20）減阻效果相對較差，Fr=0.30 時優化效果最為明顯，設計工況下（Fr=0.28）下降2.68%。

圖10 不同Fr下阻力系數對比Fig.10 Ct changing with Fr

圖11 給出了母型船和優化船在Fr為0.20～0.30 時的艏部波高對比圖。可以看出，在Fr=0.20～0.24時優化后的球鼻艏處的興波明顯低于母型船，隨著航速的增大，優化船船首處的波高與母型船的差距逐漸縮減；當Fr＞0.28 以后，優化船艏部的波形與母型船逐漸相近。優化船周圍的波形較母型船得到改善，波高低于母型船，波峰波谷影響的區域也較母型船有所減小，見圖12。

圖11 艏部波高對比Fig.11 Comparison of bow wave heights

圖12 波形圖對比Fig.12 Comparison of wave elevation contours

圖13給出了母型船與優化船在Fr為0.20～0.30時所受到壓力的對比圖?？梢钥闯觯S著Fr的增加，優化船和母型船在球鼻艏下部受力基本一致，而球鼻艏上部的差別較為明顯，優化船所受壓力相較于母型船有所降低。另外，球鼻艏首端受到的壓力隨著Fr增大而逐漸增大，優化船的球鼻艏首端壓力相對大于母型船壓力。

圖13 壓力云圖Fig.13 Pressure force contours

4 結 論

本文提出了一種利用CFD 計算與替代模型相結合的優化設計方法并成功應用于漁船球鼻艏減阻。該方法以某金槍魚船球鼻艏為研究對象，進行參數化建模，對設計變量進行靈敏度分析，尋找對總阻力系數最敏感的兩個設計變量；接著使用窮舉搜索法進行抽樣生成樣本空間，并利用CFD進行數值計算；然后通過多項式擬合生成以特征參數為自變量，總阻力系數為因變量的二元五次多項式替代模型；最后利用遺傳算法求解替代模型的最優解，獲得最優球鼻艏形狀。優化結果考慮了多個工況下的阻力性能，得出以下結論：

（1）通過靈敏度分析確定了兩個對總阻力系數最敏感的設計參數，優化后的球鼻艏形狀在垂向上變化較大，縱向上最遠端點高度較母型船降低，橫向上的變化較小；

（2）阻力分析表明，在設計工況條件下，優化后的船舶受到的興波阻力下降5.43%，粘性阻力增加了0.16%；優化后的球鼻艏在各個工況下的總阻力系數都較母型船低，總阻力平均降低2.86%，在Fr=0.22～0.30時，優化的效果較好，在設計航速（Fr=0.28）條件下，總阻力下降2.68%，Fr=0.30時優化效果最為明顯，總阻力減小3.57%；

（3）優化后船體周圍的波形較母型船得到改善，艏部波高低于母型船，波峰波谷影響的區域也較母型船有所減小。隨著Fr的增加，優化船和母型船在球鼻艏下部受力基本一致，而球鼻艏上部的差別較為明顯，優化船所受壓力相較于母型船有所降低。球鼻艏首端受到的壓力隨著Fr增大而逐漸增大，優化船的球鼻艏首端壓力相對大于母型船壓力。

通過替代模型建立球鼻艏特征參數與總阻力系數之間的關系，考慮到計算量和模型精度的影響，只對兩個對總阻力系數最敏感的設計變量進行優化。后期將會考慮更多變量之間的耦合效應，同時也會采用精度更高的替代模型，更加深入研究球鼻艏優化問題。