基于自融合方法的船體外形優(yōu)化研究

洪智超，宗 智，趙玫佳，杜金剛

（1.江蘇科技大學(xué)，江蘇鎮(zhèn)江 212100；2.大連理工大學(xué)，遼寧大連 116024；3.中國造船工程學(xué)會，北京 100861）

0 引 言

船體幾何外形優(yōu)化主要包含三個方面[1-2]：船體幾何重構(gòu)技術(shù)、船舶性能評估預(yù)報(bào)和優(yōu)化算法。設(shè)計(jì)變量受到船體幾何重構(gòu)技術(shù)的直接影響，而且設(shè)計(jì)變量是船舶性能預(yù)報(bào)和優(yōu)化的輸入?yún)?shù)，因此船體幾何重構(gòu)技術(shù)是整個優(yōu)化過程中的基礎(chǔ)。良好的船體幾何重構(gòu)技術(shù)應(yīng)該具備以下幾個優(yōu)點(diǎn)：設(shè)計(jì)變量少、幾何外形變形幅度大、生成的曲面光順以及生成的船體外形具有實(shí)用性[3]。近年來，CFD因其高精度得到了廣泛應(yīng)用，但其對計(jì)算資源的消耗也較大，因此參數(shù)數(shù)量少的船體幾何重構(gòu)技術(shù)受到了廣泛關(guān)注。Tahara[4-5]使用融合方法對一艘船的幾何外形進(jìn)行了優(yōu)化，通過對兩艘船進(jìn)行融合實(shí)現(xiàn)船體幾何重構(gòu)。該方法的優(yōu)點(diǎn)是使用的參數(shù)數(shù)量少，但船體幾何外形的變化限制其用于融合的母型船之間。特征參數(shù)方法[6]使用的參數(shù)數(shù)量也較少，與傳統(tǒng)的基于CAD 的幾何重構(gòu)技術(shù)不同；特征參數(shù)法并不使用分布在船體表面的點(diǎn)，而是使用形狀參數(shù)以及一組縱剖線來描述船體外形。該方法的缺點(diǎn)是形狀參數(shù)的選擇對設(shè)計(jì)員的經(jīng)驗(yàn)依賴性較大。

變形幅度大的幾何重構(gòu)技術(shù)由船體外形局部優(yōu)化拓展到全局優(yōu)化時，設(shè)計(jì)變量會急劇增加，從而導(dǎo)致優(yōu)化時間大幅增加，因此此類方法一般用于船體外形的局部優(yōu)化。Peri等[7]使用CFD 方法對一艘油輪進(jìn)行了優(yōu)化，使用Perturbation Surface 方法對球鼻艏進(jìn)行幾何重構(gòu)，該方法的特征是變形幅度大，但因設(shè)計(jì)變量數(shù)量的原因不適合用于全局優(yōu)化；馮梅[8]提出了一種基于徑向基函數(shù)插值（radial basis function，RBF）的曲面變形方法，并將改進(jìn)的船體曲面變形模塊嵌入團(tuán)隊(duì)開發(fā)的船舶水動力性能多學(xué)科綜合優(yōu)化平臺（SHIPMDO-WUT V6.0）；Chrismianto 和Kim[9]使用三階貝塞爾曲線、曲線平面交叉法對球鼻艏進(jìn)行了優(yōu)化；Perez等[10]使用B樣條曲面方法對球鼻艏進(jìn)行了優(yōu)化，優(yōu)化中使用CFD方法對目標(biāo)函數(shù)進(jìn)行數(shù)值預(yù)報(bào)。FFD（free form deformation）[11-12]方法也是適用于船體外形局部優(yōu)化的方法，該方法的優(yōu)點(diǎn)是幾何曲面可以任意變形。倪其軍[13]使用SBD 技術(shù)對“探索一號”科考船的艏部線型進(jìn)行優(yōu)化，其船體幾何重構(gòu)的實(shí)現(xiàn)使用FFD 方法，并利用粒子群優(yōu)化算法對線型優(yōu)化設(shè)計(jì)問題進(jìn)行求解，其總阻力減小了7.7%；馮佰威[14]使用基于徑向基插值技術(shù)的船體曲面變形方法對DTMB 5415 展開阻力優(yōu)化，興波阻力最多減小了49.66%；周廣利[15]使用NURBS 技術(shù)對DTMB 5415 船體幾何重構(gòu)展開研究，實(shí)現(xiàn)了船體高精度幾何重構(gòu)，且具有曲面變形、曲面拼接和曲面分割等功能。

本文使用自融合方法對一艘高速船展開局部及全局優(yōu)化。自融合方法由融合方法[4-5]發(fā)展而來，二者的區(qū)別在于自融合方法使用船體橫剖面，而不是將整個船體作為融合的基本元素。融合方法的基本元素為整個船體，因此需要尋找兩條以上相似船體作為融合的基本元素。而自融合方法的基本元素來源于原始船體的橫剖面，省去了尋找相似船體的麻煩。使用橫剖面作為融合基本元素會增加設(shè)計(jì)變量數(shù)量，但由此造成的數(shù)值計(jì)算量的增加有限，仍然在可接受的范圍內(nèi)。使用橫剖面作為融合基本元素能大幅增加船體外形的變化范圍，且能同時對船體的局部及全局進(jìn)行變形操作。此外，自融合方法生成的船體表面較為光滑，可直接用于后續(xù)的數(shù)值計(jì)算，無需進(jìn)行人工修正。

本文使用具有6個變量的自融合方法對高速船展開局部及全局優(yōu)化。使用CFD 方法預(yù)報(bào)船體總阻力系數(shù)（優(yōu)化的目標(biāo)函數(shù)）。CFD 方法對計(jì)算資源的消耗較大，為了降低計(jì)算量，節(jié)約計(jì)算資源，對該問題展開靈敏度分析并建立響應(yīng)面（RSM）。靈敏度分析的結(jié)果表明，全局變量（控制全船曲面變形的變量）對目標(biāo)函數(shù)的影響大幅高于局部變量（控制船體局部變形的變量）的影響。使用RSM 進(jìn)行優(yōu)化可大幅降低數(shù)值計(jì)算的工況數(shù)量，由56下降到100 以內(nèi)。本文使用MIGA（multi-island genetic algorithm）對該問題進(jìn)行優(yōu)化，與傳統(tǒng)的GA 相比，該方法能避免最優(yōu)解陷入局部陷阱，可靠性得到大幅提高。優(yōu)化后，該船的總阻力系數(shù)下降了9.76%。

1 船體自融合方法

船體自融合方法是船體幾何外形優(yōu)化的重要組成部分，該方法選取母型船的橫剖面為基本元素并使用融合因子給其賦予權(quán)重，生成新的橫剖面，最后使用新的橫剖面重構(gòu)船體三維曲面，獲得新的船型。文中使用6個自融合因子和12個由原船中選取的基礎(chǔ)橫剖面進(jìn)行船體橫剖面的融合及船體三維曲面的幾何重構(gòu)。

1.1 自融合方法基本流程

自融合方法能快速高效地修改船體幾何外形，其關(guān)鍵在于融合因子和從原始船型表面提取的橫剖面，自融合方法的計(jì)算式見式（1），基本流程見圖1。首先，根據(jù)原船基本特性及優(yōu)化目標(biāo)選取第一類橫剖面，第一類橫剖面用于后續(xù)的三維船體幾何重構(gòu)，所選取的第一類橫剖面應(yīng)盡可能完整地包含原船的外形及性能信息；然后，針對每個第一類橫剖面，從原船選取兩個第二類橫剖面，用于融合產(chǎn)生第一類橫剖面，并作為第一類橫剖面的變形空間的邊界；最后調(diào)整融合因子，產(chǎn)生新的第一類橫剖面，并進(jìn)行船體重構(gòu)獲得新船的幾何外形。自融合方法的三個重要因素有：（1）橫剖面的表達(dá)方式；（2）橫剖面的選取；（3）約束條件（本文中的約束條件為船體總排水量保持不變）。

圖1 自融合方法工作流程Fig.1 Flow chart of self-blending method

式中，CS表示橫剖面，ω表示融合因子，i、j表示橫剖面的編號。

1.2 橫剖面表達(dá)方法

橫剖面的表達(dá)應(yīng)同時具有涉及參數(shù)少和幾何外形描述精確的特性。本文通過定義橫剖面輪廓線上的節(jié)點(diǎn)，并使用樣條曲線連接各節(jié)點(diǎn)實(shí)現(xiàn)船體橫剖面的表達(dá)，其示意圖見圖2。橫剖面輪廓線上的節(jié)點(diǎn)分布方式是輪廓線構(gòu)建的關(guān)鍵，使用水線或縱剖線與輪廓線的交點(diǎn)作為節(jié)點(diǎn)是一種較為簡便的節(jié)點(diǎn)選取方式，見圖3（a）～（b），但使用這兩種方法選取節(jié)點(diǎn)會導(dǎo)致不同橫剖面上相同序號節(jié)點(diǎn)的橫坐標(biāo)或縱坐標(biāo)均為同一數(shù)值，降低了后續(xù)融合計(jì)算中節(jié)點(diǎn)移動的自由度，即降低了船體變形的幅度。本文使用輪廓線與輻射線的交點(diǎn)作為橫剖面的節(jié)點(diǎn)，見圖3（c），其融合計(jì)算式為

圖2 船體橫剖面示意圖Fig.2 Sketch of hull section

圖3 橫剖面輪廓線節(jié)點(diǎn)分布方式示意圖Fig.3 Sketch of section line nodes

式中，yN表示融合后所得到的新橫剖面的輪廓線節(jié)點(diǎn)的y坐標(biāo)，zN表示融合后所得到的新橫剖面輪廓線節(jié)點(diǎn)的z坐標(biāo)，ω表示融合因子，y表示用于融合的基礎(chǔ)橫剖面的輪廓線節(jié)點(diǎn)的y坐標(biāo)，z表示用于融合的基礎(chǔ)橫剖面的輪廓線節(jié)點(diǎn)的z坐標(biāo)，下標(biāo)i、j分別為新橫剖面、基礎(chǔ)橫剖面的序號，下標(biāo)k為橫剖面輪廓線節(jié)點(diǎn)的編號。

1.3 橫剖面的選取

基礎(chǔ)橫剖面的選取對船體外形優(yōu)化具有重要意義，直接關(guān)系到船體幾何變形的空間。從原船選取的基礎(chǔ)橫剖面應(yīng)包含船體外形的關(guān)鍵信息，尤其是與船舶性能高度相關(guān)的幾何特征。本文的研究對象為高速方尾船，其主尺度見表1，橫剖面圖見圖4，船體三維輪廓圖見圖5。該船后半段的形狀變化較小，橫剖面輪廓由兩條直線組成。因此，在0～10站之間的橫剖面中選取第0站和第9站兩個橫剖面作為第一類橫剖面。對于船體前半部分，尤其是球鼻艏部分，曲率變化較大，選取4 個橫剖面作為第一類橫剖面，全船共計(jì)6個第一類橫剖面（19.692、19.4、19、16、9、0），第二類橫剖面選取與第一類橫剖面相鄰的橫剖面。

表1 原船主尺度（模型尺度）Tab.1 Principal dimensions of original vessel(model)

圖4 船體橫剖面曲線圖Fig.4 Body plan of the vessel

圖5 船體的三維曲面圖Fig.5 3D surface of original ship hull

2 數(shù)值計(jì)算方法驗(yàn)證

本文使用CFD 方法預(yù)報(bào)船體阻力，求解器為基于有限體積法的STAR-CCM+。使用VOF 方法捕捉水和空氣的交界面，湍流模型為k-ω模型，網(wǎng)格劃分使用STAR-CCM+的切割體網(wǎng)格。此外，本文使用重疊網(wǎng)格技術(shù)模擬船體的縱搖和垂蕩運(yùn)動。數(shù)值計(jì)算的不確定度分析使用Stern 提出的收斂性分析方法[13]。阻力預(yù)報(bào)試驗(yàn)由大連理工大學(xué)完成。

本文對設(shè)計(jì)航速下（Fn=0.367）船舶阻力預(yù)報(bào)數(shù)值模擬的迭代收斂性和網(wǎng)格收斂性展開分析。迭代收斂性通過總阻力系數(shù)的時歷曲線進(jìn)行分析，網(wǎng)格收斂性分析以三套不同尺度網(wǎng)格的數(shù)值計(jì)算結(jié)果為基礎(chǔ)?？傋枇ο禂?shù)的時歷曲線見圖6。迭代不確定度為時歷曲線最后幾個周期的振幅，即為時歷曲線收斂后周期性振蕩的幅值，為0.28%。不同尺度網(wǎng)格的數(shù)量及阻力預(yù)報(bào)結(jié)果見表2，網(wǎng)格尺寸的變化率為 2，該值的選取參照Stern[16]的論文。網(wǎng)格收斂性驗(yàn)證的結(jié)果見表3。網(wǎng)格收斂因子RG為0.123，小于1，即網(wǎng)格收斂性符合單調(diào)收斂。誤差E為1.45%S，網(wǎng)格不確定度UG為0.37%，迭代不確定度為0.28%，模擬不確定度為0.22%，均為較低水平，說明使用該方法進(jìn)行阻力預(yù)報(bào)是可行的。

圖6 原始船型的CT時歷曲線Fig.6 Time history of CT of original vessel

表2 不同網(wǎng)格尺度下原始船型的總阻力系數(shù)（Fn=0.367）Tab.2 Results of CTunder different grid sizes（Fn=0.367）

表3 網(wǎng)格收斂性分析（S為網(wǎng)格1尺度下總阻力系數(shù)CT的數(shù)值計(jì)算結(jié)果）Tab.3 Grid convergence analysis(S refers to the simulation result of CT under Grid 1)

3 CFD計(jì)算結(jié)果

本文開展了一系列數(shù)值計(jì)算以構(gòu)建關(guān)聯(lián)設(shè)計(jì)變量和CT的響應(yīng)面關(guān)系式。在討論響應(yīng)面之前首先介紹一些基本概念。數(shù)值模擬空間是一個由設(shè)計(jì)變量的最小及最大值定義的范圍，變量的水平是指變量的不同取值，響應(yīng)是指總阻力系數(shù)CT的值。使用響應(yīng)面進(jìn)行優(yōu)化的流程見圖7。首先，使用靈敏度分析研究各設(shè)計(jì)變量對響應(yīng)的影響，并以此為依據(jù)定義響應(yīng)面的數(shù)值模擬空間。其次，使用中心復(fù)合設(shè)計(jì)（CCD）方法確定數(shù)值模擬計(jì)算工況并展開計(jì)算，根據(jù)結(jié)果構(gòu)建響應(yīng)面，并評估其可靠性。最后，以響應(yīng)面為基礎(chǔ)展開優(yōu)化，獲得最優(yōu)解。

圖7 使用響應(yīng)面的流程圖Fig.7 Flow chart of the application of RS

3.1 靈敏度分析結(jié)果

本文使用L25（56）正交試驗(yàn)表及直觀分析方法對設(shè)計(jì)變量展開靈敏度分析，共包含25 個計(jì)算工況、6 個因子（設(shè)計(jì)變量），每個因子包含5個水平（均布在設(shè)計(jì)空間）。6個因子及靈敏度分析的數(shù)值模擬空間見表4，其中的4個設(shè)計(jì)變量（ω3，ω4，ω5，ω6）及其對應(yīng)的橫剖面19、16、9、0 用于控制船體的全局形狀變化，3 個設(shè)計(jì)變量（ω1，ω2，ω3）及其對應(yīng)的橫剖面19.672、19.4、19用于船體的局部外形（球鼻艏）控制。

表4 設(shè)計(jì)變量空間Tab.4 Simulation domain of the design variables

正交試驗(yàn)表中所有工況的數(shù)值計(jì)算結(jié)果見圖8，由圖中的結(jié)果可知，25個工況中CT的波動幅度為22.98%。對正交試驗(yàn)的結(jié)果進(jìn)行直觀分析，其結(jié)果見圖9，圖中的結(jié)果表明因子ω6對總阻力系數(shù)的影響最大，說明方尾的形狀對總阻力系數(shù)的影響最為顯著。因子ω1、ω2、ω3對總阻力系數(shù)的影響遠(yuǎn)小于因子ω4、ω5、ω6的影響，說明局部形狀的變化對總阻力系數(shù)的影響小于全局形狀變化的影響，這與實(shí)際情況相符。根據(jù)圖9的結(jié)果，最優(yōu)解在融合因子的水平為[1，4，2，3，5，5]的鄰域內(nèi)。

圖8 正交試驗(yàn)數(shù)值計(jì)算結(jié)果Fig.8 Calculation results of the orthogonal test

圖9 設(shè)計(jì)變量直觀分析圖Fig.9 Main effect of screening design

3.2 響應(yīng)面結(jié)果

直觀分析能給出各因子對目標(biāo)函數(shù)的影響趨勢，但無法考慮因子間的相互影響，因此本文以靈敏度分析為基礎(chǔ)，在其得出的最優(yōu)解的鄰域內(nèi)進(jìn)行CCD數(shù)值試驗(yàn)設(shè)計(jì)，CCD 數(shù)值試驗(yàn)設(shè)計(jì)中各因子的取值范圍見表4，每個因子取5 個水平，共計(jì)77 個工況，其數(shù)值計(jì)算結(jié)果見圖10。根據(jù)CCD 的數(shù)值計(jì)算結(jié)果構(gòu)建一個基于二次多項(xiàng)式的響應(yīng)面，該響應(yīng)面建立了設(shè)計(jì)變量與總阻力系數(shù)之間的映射關(guān)系，且該響應(yīng)面包含了因子間的相互影響，響應(yīng)面的計(jì)算式為

圖10 CCD數(shù)值計(jì)算結(jié)果Fig.10 Results of CT by CCD

為檢驗(yàn)響應(yīng)面的有效性，檢驗(yàn)其是否能真實(shí)反映設(shè)計(jì)變量與目標(biāo)函數(shù)之間的關(guān)系，對響應(yīng)面展開方差分析，其結(jié)果見表5（表中SS為平方和，DF為自由度，MS為均方）。方差分析主要用于評估響應(yīng)面，檢驗(yàn)其是否能真實(shí)反映設(shè)計(jì)變量與目標(biāo)函數(shù)之間的關(guān)系。方差分析的中心思想是通過對比由設(shè)計(jì)變量變化引起的目標(biāo)函數(shù)波動與由隨機(jī)誤差引起的目標(biāo)函數(shù)波動來判斷響應(yīng)面的結(jié)果是否顯著。文中響應(yīng)面的F值為21.9，說明該響應(yīng)面是顯著的，即該響應(yīng)面是有效的。

表5 響應(yīng)面的方差分析Tab.5 Variance analysis of RS

4 船體總阻力系數(shù)優(yōu)化

本文以控制船體曲面變形的融合因子為設(shè)計(jì)變量對船舶的總阻力系數(shù)CT進(jìn)行優(yōu)化，其數(shù)學(xué)描述見式（4）。總阻力系數(shù)由式（3）定義，該式為以CCD 的數(shù)值計(jì)算結(jié)果為基礎(chǔ)構(gòu)建的響應(yīng)面，其設(shè)計(jì)變量為6個融合因子。響應(yīng)面在水平2 到水平4 之間具有較高的精度，因此本文的優(yōu)化空間即限定在水平2和水平4之間的設(shè)計(jì)空間內(nèi)，各因子的具體范圍見表4。優(yōu)化的約束條件是保持船體排水體積不變。優(yōu)化算法為多島遺傳算法，多島遺傳算法的各參數(shù)見表6。初始種群的總阻力系數(shù)均布在0.005 441到0.005 594之間。

表6 MIGA參數(shù)Tab.6 MIGA parameters

式中，gi為約束條件。

多島遺傳算法具有一定的偶然性，本文對該優(yōu)化問題進(jìn)行10 次優(yōu)化計(jì)算，結(jié)果表明10 次優(yōu)化的最優(yōu)解極為接近，設(shè)計(jì)變量之間的誤差小于0.05%，目標(biāo)函數(shù)之間的誤差小于0.000 1%。其中一個優(yōu)化計(jì)算的結(jié)果見圖11 和圖12。結(jié)果表明，經(jīng)過20 代的進(jìn)化之后，總阻力系數(shù)最終收斂到0.005 427。最優(yōu)解所對應(yīng)的各設(shè)計(jì)變量值見表7。為驗(yàn)證該最優(yōu)解，本文使用CFD 方法對該最優(yōu)解進(jìn)行數(shù)值計(jì)算，總阻力系數(shù)CT為0.005 422，與基于響應(yīng)面的最優(yōu)解吻合良好。優(yōu)化后，總阻力系數(shù)下降9.76%。最優(yōu)船型與原始船型的對比見圖13，由圖中的結(jié)果可知，優(yōu)化后型寬、型深及球鼻艏變化較小，船體折邊線垂向位置升高且船長略有增加。原船與優(yōu)化后船型的波形對比見圖14，由圖可知，優(yōu)化后船艏興波小于原船。

圖11 CT迭代過程Fig.11 Changes of CT in optimization

圖12 CT迭代曲線Fig.12 Convergent results of CT

圖13 最優(yōu)船型和原始船型的對比圖Fig.13 Comparison of geometies between original and optimized ship hull forms

圖14 原船與最優(yōu)船波形圖對比（上半部為原船，下半部分為最優(yōu)船）Fig.14 Wave patterns of original and optimized vessels(upper half:original hull;lower half:optimized hull)

表7 最優(yōu)解Tab.7 Optimized solution

5 結(jié) 語

本文提出了一種船體變形方法——自融合方法，其核心思想是從原始船型中提取特征橫剖面用于融合操作，生成新的橫剖面并進(jìn)行三維船體的幾何重構(gòu)，產(chǎn)生新的船型，具有使用變量少、生成的船體外形光順等優(yōu)點(diǎn)。使用6個融合因子和12個從原船提取的基準(zhǔn)橫剖面，結(jié)合CFD方法、響應(yīng)面方法和多島遺傳算法，以總阻力系數(shù)最小為優(yōu)化目標(biāo)，排水量保持不變?yōu)榧s束條件，對一艘高速方尾船展開局部及全局優(yōu)化，即對其全船外形和球鼻艏形狀展開優(yōu)化。優(yōu)化后該船在保持排水量不變的前提下，設(shè)計(jì)航速狀態(tài)下的總阻力系數(shù)下降了9.76%，證明了該方法的工程應(yīng)用價(jià)值。