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

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

0 引 言

近年來，由于過度捕撈及環(huán)境污染導(dǎo)致近海漁業(yè)資源匱乏。為緩解我國漁業(yè)壓力、建設(shè)新型現(xiàn)代化漁業(yè)體系、推進海洋強國戰(zhàn)略部署，發(fā)展遠(yuǎn)洋漁業(yè)勢在必行。2020年“中國遠(yuǎn)洋漁業(yè)履約白皮書”提出支持漁船更新改造，鼓勵采用安全高效、節(jié)能減排的技術(shù)和船用裝備。船舶參數(shù)化水動力構(gòu)型為綠色節(jié)能船型的研發(fā)提供了極大的便利，是船舶行業(yè)實現(xiàn)“碳達峰”、“碳中和”目標(biāo)的技術(shù)途徑之一。球鼻艏可在一定程度上減小船舶阻力，但其泛化性較弱，即同一球鼻艏應(yīng)用到其他船時并不一定能夠獲得滿意的減阻效果。因此，球鼻艏形狀、參數(shù)等與船體有恰當(dāng)?shù)呐浜鲜乔虮囚挤桨冈O(shè)計、優(yōu)化的重點方向。

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

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

1 基本理論與方法

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

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

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

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

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

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

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

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

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

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

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

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

2 問題描述

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

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

圖2 網(wǎng)格劃分Fig.2 Mesh division

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

圖3 KCS阻力數(shù)據(jù)對比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，主要船型參數(shù)如表1所示。在球鼻艏建模過程中，針對球鼻艏主要參數(shù)及幾何特征，共有7個球鼻艏特征參數(shù)，其初始值及取值范圍如表2所示。表2中，X1和X2是與Abl相關(guān)的特征參數(shù)，X4與Abt相關(guān)。同時還設(shè)置了特征曲線，它們分別為最大寬度曲線以及最大寬度高度曲線，即Bb在船長方向上，從最大橫剖面處至球鼻艏最遠(yuǎn)點寬度和高度的變化曲線。

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

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

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

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

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

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

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

其中擬合系數(shù)如表3所示。

表3 擬合系數(shù)Tab.3 Fitting coefficients

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

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

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

引入替代模型后，優(yōu)化問題可表示為

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

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

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

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

3 結(jié)果與討論

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

表4 優(yōu)化結(jié)果Tab.4 Optimization results

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

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

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

圖10 為母型船與優(yōu)化船在Fr為0.20～0.30 工況下的總阻力系數(shù)對比。圖中顯示，在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%。可以看出，本文優(yōu)化后的球鼻艏在低航速時（Fr≤0.20）減阻效果相對較差，F(xiàn)r=0.30 時優(yōu)化效果最為明顯，設(shè)計工況下（Fr=0.28）下降2.68%。

圖10 不同F(xiàn)r下阻力系數(shù)對比Fig.10 Ct changing with Fr

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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