改進(jìn)的徑向基插值方法在船型優(yōu)化中的應(yīng)用

馮佰威 王首茗 馮梅

(1.武漢理工大學(xué) 高性能船舶技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430063； 2.武漢理工大學(xué) 船海與能源動力工程學(xué)院，湖北 武漢 430063； 3.國家知識產(chǎn)權(quán)局專利局專利審查協(xié)作湖北中心，湖北 武漢 430063)

傳統(tǒng)的船舶設(shè)計(jì)大多采用母型船改造法，該方法經(jīng)常得到的是可行方案，而非最優(yōu)方案。船型優(yōu)化將計(jì)算流體力學(xué)與最優(yōu)化理論結(jié)合并運(yùn)用在船型設(shè)計(jì)過程中，可獲得在一定約束條件下的性能最優(yōu)的船型。

船體曲面修改是船型優(yōu)化的基礎(chǔ)。眾多學(xué)者對船體曲面變形開展了深入研究，目前常用變形方法有自由變形方法、融合方法、徑向基函數(shù)(RBF)變形方法以及參數(shù)化方法。李冬琴等[1]將非均勻有理B樣條(NURBS)基函數(shù)和非均勻控制點(diǎn)框架結(jié)合，建立適用于船體幾何的基于NURBS基函數(shù)的自由變形技術(shù)(NFFD)，以某CNG運(yùn)輸船為例完成球艏、船首和船尾部分的自動變形。馮佰威[2]在船體曲面NURBS表達(dá)的基礎(chǔ)上，使用融合方法對DTMB5415聲吶罩區(qū)域和1300TEU集裝箱船球鼻艏區(qū)域進(jìn)行了興波阻力優(yōu)化。Tahara等[3]使用融合方法將兩種尾部型線差別較大的船型融合成新船型，避免了不符合實(shí)際的船型的出現(xiàn)。Kim等[4]將修改橫剖面面積曲線方法和徑向基函數(shù)插值方法相結(jié)合，達(dá)到對KCS船的船體曲面進(jìn)行全局和局部修改的目的；沈通等[5]采用NURBS表達(dá)與RBF相結(jié)合的方法，實(shí)現(xiàn)了船體曲面的三維自動變形，研究結(jié)果表明，基于RBF的船體曲面修改方法適用于船型優(yōu)化，符合工程實(shí)際需要；周昊等[6]采用RBF方法完成了DTMB5415船型的艏部型線優(yōu)化。馮佰威等[7]采用RBF方法，完成了不同弗勞德數(shù)下的船舶阻力性能優(yōu)化。張萍[8]提出了常規(guī)船型的參數(shù)化設(shè)計(jì)方法，由特征參數(shù)設(shè)計(jì)特征曲線，生成橫剖面曲線，得到光順的船體曲面。董楓[9]建立有效的船體參數(shù)化模型，針對復(fù)雜船型的各類橫剖線構(gòu)造了完整的形狀參數(shù)集，采用改進(jìn)的曲面放樣法來構(gòu)造船體曲面模型。胡春平[10]用Friendship軟件對KCS船型進(jìn)行全參數(shù)化建模，應(yīng)用Shipflow軟件以興波阻力為目標(biāo)對型線進(jìn)行了優(yōu)化。

文獻(xiàn)[5]中對各種船體曲面變形方法進(jìn)行了詳細(xì)的對比分析。在諸多的船體曲面變形方法中，由于RBF船體曲面變形方法具有靈活高效，變形參數(shù)可控的特點(diǎn)，因此具有較強(qiáng)的工程應(yīng)用價(jià)值。

在基于RBF的船體曲面變形方法中，支撐半徑為關(guān)鍵參數(shù)，影響了船體曲面變形范圍的大小。楊烜等[11]研究了徑向基函數(shù)緊支撐集選取的問題，對任意標(biāo)志點(diǎn)集，采用構(gòu)造Delaunay三角剖分方法來確定基函數(shù)支撐集的大??；Zhang等[12]使用具有多個支撐半徑的緊支撐徑向基函數(shù)應(yīng)用于心腦圖像的配準(zhǔn)中。

本團(tuán)隊(duì)在前期研究中采用Delaunay三角剖分方法求取支撐半徑，并成功應(yīng)用于船體曲面的光順變形。但尚存在的問題是：在某些特殊情況下，Delaunay三角剖分無法獲得支撐半徑，限制了該方法的進(jìn)一步推廣應(yīng)用。為解決上述問題，文中在前期研究基礎(chǔ)上加以改進(jìn)，提出了支撐半徑的動態(tài)求取方法，并將該方法應(yīng)用到標(biāo)準(zhǔn)船型優(yōu)化中，以驗(yàn)證其可行性。

1 基于徑向基函數(shù)插值的曲面變形 方法

1.1 徑向基函數(shù)及其插值技術(shù)

徑向基函數(shù)是一種沿徑向?qū)ΨQ的標(biāo)量函數(shù)，以空間中任一點(diǎn)Xi到某一中心x之間的歐氏距離r=‖X-Xi‖為自變量。各基函數(shù)的具體形式為[13]：

φ(‖X-Xi‖),i=1,2,…,m

(1)

在船體曲面變形中應(yīng)用RBF方法時(shí)，使用完整形式的插值方程：

(2)

船體曲面上中心點(diǎn)X=(x,y,z)移動的距離用s(X)表示，點(diǎn)的數(shù)量用n表示，Xi為任一點(diǎn)，p(X)為低階多項(xiàng)式。其形式為

p(X)=c1x+c2y+c3z+c4

(3)

得到完整插值方程的關(guān)鍵在于求得λ和c。該方法的詳細(xì)介紹請參考文獻(xiàn)[13]。

1.2 支撐半徑求取存在的問題

在船體曲面變形中，選取WendlandΨ3,1作為徑向基函數(shù)插值的基函數(shù)，基函數(shù)選擇參考文獻(xiàn)[13]中的情況，該基函數(shù)形式如下：

(4)

在該函數(shù)中添加支撐半徑R，使基函數(shù)變?yōu)槿缦滦问剑?/p>

φ(r)=(1-r/R)4(4r/R+1)

(5)

支撐半徑對基函數(shù)的影響范圍進(jìn)行控制，決定了曲面變形的質(zhì)量。本團(tuán)隊(duì)主要采用Delaunay三角剖分的方法實(shí)現(xiàn)支撐半徑的求取[14]。該方法尚存在的問題是：只考慮了可以進(jìn)行三角剖分的情況，對于無法進(jìn)行三角剖分的情況未作說明，且運(yùn)用三角剖分求解支撐半徑僅考慮到點(diǎn)的分布對支撐半徑的影響，點(diǎn)的坐標(biāo)變化量對支撐半徑的影響并未考慮。這些缺點(diǎn)限制了船體曲面變形的應(yīng)用范圍。

2 支撐半徑的動態(tài)求取

2.1 支撐半徑求取分類

支撐半徑動態(tài)求取的思路是：以可變點(diǎn)的數(shù)量為依據(jù)，要使變形后的船體曲面不出現(xiàn)鞍點(diǎn)且保持拓?fù)?，需要綜合考慮可變點(diǎn)的分布和可變點(diǎn)的坐標(biāo)變化量對支撐半徑的影響。在采用Delaunay三角剖分求取支撐半徑時(shí)，重點(diǎn)在于可變點(diǎn)是否能構(gòu)成三角形，著重考慮的是可變點(diǎn)的分布對支撐半徑的影響?？勺凕c(diǎn)的分布情況可分為圖1中所示的幾類，對于n<3或者n≥3且所有點(diǎn)共線，此時(shí)僅考慮可變點(diǎn)的分布情況，考慮無法采用Delaunay三角剖分時(shí)支撐半徑的求取。同時(shí)將可變點(diǎn)坐標(biāo)變化量考慮在內(nèi)。

圖1 可變點(diǎn)分布情況

2.1.1 考慮可變點(diǎn)坐標(biāo)變化量的支撐半徑求取

在船體曲面變形中要求保持拓?fù)?。必要條件是函數(shù)S是連續(xù)的，并且雅可比矩陣的行列式必須在圖像的每個點(diǎn)處為正：

det(?s)>0

(6)

式中，?s為雅可比矩陣。

為了證明變換的全局不變性屬性，使用如下的定理[15]。

定理1 如果s:Ω→Rn，其中Ω是Rn中一個封閉的矩形區(qū)域，是可微分映射，使得雅可比矩陣?s是Ω中所有x的P矩陣，那么s在Ω中是唯一的。

注意，如果n×n實(shí)矩陣的所有主子式都是正的，則稱其為P矩陣。

為了變換函數(shù)s，分析了單個點(diǎn)p的雅可比矩陣，它包括一個恒等映射和一個局部基函數(shù)的疊加。對于2維的單個點(diǎn)有：

s1(X)=x+Δ1φ(‖X-P‖)

s2(X)=y+Δ2φ(‖X-P‖)

(7)

式中，φ是局部RBF基函數(shù)，Δ1和Δ2分別是在x和y坐標(biāo)方向上從初始點(diǎn)p到目標(biāo)點(diǎn)q的位移，并且p是任意位置的點(diǎn)。將式(7)代入式(6)，所得如下：

式中，φ表示φ(‖X-P‖)，r=‖X-P‖，如果令Δ=max(Δ1,Δ2)，在γ=π/4時(shí)，得到二維下的最壞情況：

(8)

對于三維的情況，獲得了類似的條件：

如果再次設(shè)置Δ=max(Δ1,Δ2,Δ3)，獲得在三維中最差的情況：

(9)

利用式(8)和式(9)，可以顯示雅可比行列式的所有主子式分別對于二維和三維情況是正的。由此得出，如果式(8)和式(9)成立，則變換保留拓?fù)洹?φ/?r的最小值取決于支撐半徑R。表1總結(jié)了不同局部基函數(shù)的?φ/?r。在表1和式(8)、式(9)的前提下，可以導(dǎo)出R的條件。如表2所示。

表1 ?φ/?r的最小值

表2 給定位于Δ時(shí)最小的R

這些條件僅適用于在支撐半徑R內(nèi)沒有其他點(diǎn)的孤立點(diǎn)。對于具有相交支撐區(qū)域的點(diǎn)，R的最小值還取決于點(diǎn)p的位置，這導(dǎo)致更復(fù)雜的計(jì)算。表2為選擇支撐半徑R的求取提供了線索。

在單標(biāo)志點(diǎn)情況下，為了保持變換后的拓?fù)潢P(guān)系，WendlandΨ3,1基函數(shù)的支撐半徑在三維時(shí)應(yīng)滿足R>3.66Δ，Δ=max(Δ1,Δ2,Δ3)，即支撐半徑與點(diǎn)的坐標(biāo)變化量有關(guān)，說明支撐半徑有下界。

2.1.2 無法進(jìn)行三角剖分時(shí)的支撐半徑求取

無法進(jìn)行三角剖分的情況如下：

(1)可變點(diǎn)的數(shù)量小于3。Delaunay三角剖分的基本思想是以三角網(wǎng)格的形式表征某一點(diǎn)與周圍鄰近點(diǎn)的拓?fù)潢P(guān)系，當(dāng)可變點(diǎn)數(shù)量小于3個時(shí)無法構(gòu)成三角形。

(2)可變點(diǎn)數(shù)量不小于3，但這些點(diǎn)共線?；A(chǔ)三角網(wǎng)增長法在尋找點(diǎn)和擴(kuò)展邊的時(shí)候要滿足相應(yīng)的準(zhǔn)則，當(dāng)不滿足準(zhǔn)則時(shí)無法進(jìn)行三角剖分。

針對上述情況，給出的解決方法如下：

(1)可變點(diǎn)的數(shù)量為1時(shí)，僅考慮可變點(diǎn)坐標(biāo)變化量的影響，選取的支撐半徑為R=3.66Δ。

(2)可變點(diǎn)的數(shù)量不小于2時(shí)，此時(shí)為所有點(diǎn)共線的情況。對于WendlandΨ3,1基函數(shù)，楊烜等[11]已經(jīng)證明：在二維兩個控制點(diǎn)的情況下，只要保證兩控制點(diǎn)間的區(qū)域不出現(xiàn)鞍點(diǎn)，就可以保證在其它位置不出現(xiàn)鞍點(diǎn)。并且當(dāng)支撐半徑取兩點(diǎn)距離的2倍時(shí)，兩點(diǎn)連線中點(diǎn)處的二階導(dǎo)數(shù)小于0，不出現(xiàn)鞍點(diǎn)。記相鄰兩點(diǎn)間最大的歐氏距離為rd，僅考慮可變點(diǎn)分布的影響，取支撐半徑為R=2rd。綜合考慮可變點(diǎn)的分布和坐標(biāo)變化量的影響，最終取支撐半徑為R=max{2rd,3.66Δ}。

2.2 支撐半徑的動態(tài)求取框架

支撐半徑動態(tài)求取的具體算法框架如圖2所示，詳細(xì)的步驟可總結(jié)如下：

(1)首先在特定文件中讀取出可變點(diǎn)的個數(shù)n(n為正整數(shù))，然后求出各個點(diǎn)在每個方向(x,y,z)上的坐標(biāo)變化量，并且求取出各變化量的最大值Δ。

(2)判斷一 用于確定可變點(diǎn)的數(shù)量n是否小于3。若n小于3，則因?yàn)辄c(diǎn)的個數(shù)太少無法構(gòu)成三角形，從而不能進(jìn)行Delaunay三角剖分，并且需要進(jìn)行判斷二；若n不小于3，則需要做進(jìn)一步的判斷，即判斷三。

(3)判斷二 用于將可變點(diǎn)個數(shù)小于3的情況分為兩種：當(dāng)n等于1時(shí)，不用考慮可變點(diǎn)的分布對支撐半徑的影響，支撐半徑取為支撐半徑的下界，即R=3.66Δ；當(dāng)n等于2時(shí)，考慮可變點(diǎn)分布的影響，從而求取相鄰兩點(diǎn)間的最大距離rd，之后需要進(jìn)行判斷四。

(4)判斷三 用于將可變點(diǎn)個數(shù)不小于3的情況分為兩種：當(dāng)所有可變點(diǎn)共線的時(shí)候，不滿足三角形的構(gòu)成條件，無法進(jìn)行Delaunay三角剖分，考慮可變點(diǎn)分布的影響，求取相鄰可變點(diǎn)間的最大距離rd，之后進(jìn)行判斷四；當(dāng)不是所有可變點(diǎn)都共線的時(shí)候，考慮可變點(diǎn)分布的影響，采用Delaunay三角剖分求取三角形的最大邊長rd，之后進(jìn)行判斷四。

(5)判斷四 用于對可變點(diǎn)的分布位置和可變點(diǎn)的坐標(biāo)變化量進(jìn)行綜合考量，同時(shí)滿足變形之后不能產(chǎn)生鞍點(diǎn)，且保持拓?fù)浣Y(jié)構(gòu)，最終選取2rd和3.66Δ中較大的為支撐半徑R。

圖2 支撐半徑動態(tài)算法框架

2.3 支撐半徑的動態(tài)求取實(shí)例

將支撐半徑的動態(tài)求取用于具體船型的變化，首先考慮可變點(diǎn)數(shù)量為2時(shí)的情況，該情況與可變點(diǎn)三點(diǎn)共線時(shí)求取策略相同。在DTMB5415船型的艏部選擇兩個可變點(diǎn)Y1、Y2，兩點(diǎn)在設(shè)計(jì)水線附近，且坐標(biāo)均在y軸方向變化，具體坐標(biāo)數(shù)據(jù)見表3，經(jīng)計(jì)算其3.66Δ為0.039 2，2rd為0.915 6，最終支撐半徑為0.915 6；選擇甲板邊線以及聲吶球艏處的點(diǎn)為固定點(diǎn)，如圖3所示。運(yùn)用支撐半徑動態(tài)求取算法求取相應(yīng)的支撐半徑R，并且采用徑向基函數(shù)插值方法得到變形后的船型，初始船和變形船的型線對比如圖4所示?？梢钥闯鲎冃未?水線(吃水為0處的水線)以下的型線未發(fā)生變化；在0水線以上的部分由橫剖線圖可知，其略微變肥，且首部型線內(nèi)凹；半寬水線(首部)顯示甲板邊線沒有發(fā)生變化，部分水線變凹。

表3 船型算例坐標(biāo)數(shù)據(jù)表

圖3 可變點(diǎn)數(shù)量為2時(shí)選點(diǎn)位置圖

圖4 可變點(diǎn)數(shù)量為2時(shí)型線對比圖

考慮可變點(diǎn)數(shù)量為1時(shí)，在S60船型艏部選取一個可變點(diǎn)，可變點(diǎn)沿Y軸方向變化?？勺凕c(diǎn)和固定點(diǎn)位置如圖5所示，具體坐標(biāo)數(shù)據(jù)如表4所示。運(yùn)用支撐半徑動態(tài)求取算法求取相應(yīng)的支撐半徑，由于可變點(diǎn)數(shù)量為1，因此支撐半徑值為3.66Δ，經(jīng)計(jì)算3.66Δ為0.023，因此其支撐半徑為0.023。并且采用徑向基函數(shù)插值方法得到變形后的船型，初始船和變形船的型線對比如圖6所示，首部橫剖線形狀外凸，部分縱剖線呈S狀。

圖5 可變點(diǎn)數(shù)量為1時(shí)選點(diǎn)位置圖

表4 船型算例坐標(biāo)數(shù)據(jù)表

圖6 可變點(diǎn)數(shù)量為1時(shí)型線對比圖

綜上所述，支撐半徑動態(tài)求取方法可以進(jìn)行特殊情況下支撐半徑的求取，且所求取的支撐半徑可以應(yīng)用于徑向基函數(shù)插值，實(shí)現(xiàn)船體曲面變形。

3 應(yīng)用實(shí)例

3.1 優(yōu)化問題描述

將改進(jìn)的徑向基插值曲面變形方法應(yīng)用于SHIPMDO-WUT平臺，對S60船型進(jìn)行曲面的變形優(yōu)化。用于變形優(yōu)化的S60船模主尺度信息如表5所示。

表5 船模主尺度信息

優(yōu)化目標(biāo)定義如下(Fn為傅汝德數(shù)):

fobj=minCw,Fn=0.30

(10)

其中興波阻力系數(shù)計(jì)算公式為

(11)

式中，ρ為海水的密度，V為航速，S為濕表面積，Rw為興波阻力。

1)變量選擇

選取如圖7所示的4個曲面上的點(diǎn)，且每個點(diǎn)在y、z兩個方向變化，例如點(diǎn)1的y坐標(biāo)為設(shè)計(jì)變量Y1，點(diǎn)1的z坐標(biāo)為設(shè)計(jì)變量Z1，依次類推，共4個可變點(diǎn)，8個坐標(biāo)變化，即8個設(shè)計(jì)變量。具體設(shè)計(jì)變量及變量范圍如表6所示。

圖7 選點(diǎn)位置圖

表6 設(shè)計(jì)變量及其范圍

2)約束條件

靜水力約束：優(yōu)化船排水體積不得小于母型船。

型線約束：各主尺度要素船長、船寬、吃水均保持不變；首部輪廓線，基線保持不變，船體部分站位保持不變。固定點(diǎn)位置如圖7所示。

3)優(yōu)化算法

選擇多目標(biāo)粒子群算法(MOPSO)修改優(yōu)化變量，設(shè)置20代，每代粒子數(shù)為50。在粒子群優(yōu)化算法中，粒子速度和粒子位置更新見下式[16]。

(12)

式中:vi(t+1)為第t+1代、第i個粒子的速度;xi(t+1)為第t+1代、第i個粒子的位置；w為權(quán)重因子；c1、c2為學(xué)習(xí)因子；r1、r2表示0到1的隨機(jī)數(shù);α表示局部最優(yōu)解;β表示全局最優(yōu)解。當(dāng)不滿足優(yōu)化收斂準(zhǔn)則時(shí)，根據(jù)上式更新優(yōu)化變量。

4)參數(shù)化變形模塊

本流程的參數(shù)化變形模塊即為基于徑向基函數(shù)插值的船體曲面變形程序。該模塊在迭代過程中讀取更新的優(yōu)化變量，通過徑向基函數(shù)的插值操作，可實(shí)現(xiàn)船體曲面的光順變形。

具體的優(yōu)化流程如圖8所示。

圖8 優(yōu)化流程圖

3.2 優(yōu)化結(jié)果分析

優(yōu)化所得的1050組結(jié)果如圖9所示。圖中標(biāo)注了初始船的興波阻力系數(shù)位置，選擇其中興波阻力系數(shù)最小的船型為優(yōu)化船。表7、表8分別列出了初始船和優(yōu)化船的靜水力數(shù)據(jù)、水動力數(shù)據(jù)。為了節(jié)省時(shí)間，優(yōu)化過程中采用Shipflow軟件計(jì)算船體興波阻力系數(shù)；為了驗(yàn)證優(yōu)化效果，優(yōu)化結(jié)束后采用Fine/Marine軟件分別計(jì)算初始船和優(yōu)化船的總阻力。優(yōu)化船相較初始船而言，排水體積增加0.83%，浮心縱向位置減小0.38%，濕表面積增加0.62%，興波阻力減小10.46%，總阻力減小3.24%。

圖9 優(yōu)化收斂圖

表7 靜水力數(shù)據(jù)

表8 水動力數(shù)據(jù)

由式(11)可知，在排水體積和浮心縱向位置均滿足約束的條件下，且航速一定時(shí)，優(yōu)化船的濕表面積增加，興波阻力系數(shù)減小，說明興波阻力減小，且其減小是由船體型線的變化導(dǎo)致的；優(yōu)化船總阻力減小，說明在優(yōu)化時(shí)達(dá)到了通過興波阻力系數(shù)減小從而使總阻力減小的目的。

圖10顯示的是初始船和優(yōu)化船的型線對比圖，圖11為初始船和優(yōu)化船的三維模型圖。波形圖和波切圖的對比分別在圖12、圖13中顯示。圖13中LY為Y方向波切位置距船體中縫剖面的距離，圖14為壓力分布圖。從型線對比圖中明顯可以看到型線變化主要在首部，產(chǎn)生了隱形球艏，部分縱剖面呈S型，對于高速時(shí)產(chǎn)生有利的興波干擾有比較大的作用。由波切圖可以看出，波高在波峰和波谷位置處的幅值減少，因此船首到船中部分(x/L為0～0.5，L為船長)的興波幅值略微減小，船中到船尾部分(x/L為0.5～1.0)的興波幅值明顯減小，船尾及尾后0.5L部分(x/L為1.0～1.5)興波幅值略微減小。說明隱形球艏的產(chǎn)生對船體興波產(chǎn)生了有利干擾，使得興波阻力減小。從壓力分布圖可以看出，優(yōu)化船船首部分較初始船壓力增加，船尾部分變化較小，因此其首尾壓差阻力增加。這是由于船首產(chǎn)生了隱形球艏，濕表面積增加，黏性阻力增大。由于S60船型為中高速船，興波阻力在總阻力中占比較大，而黏性阻力占比較小。因此優(yōu)化船興波阻力有明顯降低，總阻力降低。

圖10 優(yōu)化前后型線對比圖

圖11 三維模型圖

圖12 波形圖

圖13 波切圖

圖14 壓力分布圖

4 結(jié)論

支撐半徑是徑向基函數(shù)插值中的重要參數(shù)，對船體曲面變形有較大影響。本文重點(diǎn)針對支撐半徑的求取進(jìn)行了改進(jìn)。并成功應(yīng)用于S60船型的優(yōu)化方案中。結(jié)論如下：

(1)改進(jìn)后的支撐半徑求取既考慮了可變點(diǎn)位置分布的影響，又考慮了可變點(diǎn)坐標(biāo)變化量的影響，實(shí)現(xiàn)了支撐半徑的動態(tài)求取。

(2)可變點(diǎn)的數(shù)量為1時(shí)，僅考慮可變點(diǎn)坐標(biāo)變化量的影響，選取的支撐半徑為R=3.66Δ。

(3)可變點(diǎn)的數(shù)量不小于2時(shí)，此時(shí)為所有點(diǎn)共線的情況。取支撐半徑為R=max{2rd,3.66Δ}。

(4)可變點(diǎn)數(shù)量不小于3且不共線時(shí)，采用Delaunay三角剖分求取三角形的最大邊長rd，最終選取2rd和3.66Δ中較大的為支撐半徑R。

綜上所述，動態(tài)支撐半徑的求取方法運(yùn)用到基于徑向基函數(shù)的船型優(yōu)化中是可行的，具有一定的工程應(yīng)用價(jià)值。