改进的径向基插值方法在船型优化中的应用

2022-04-29冯佰威王首茗冯梅

华南理工大学学报（自然科学版） 2022年3期

冯佰威王首茗冯梅

(1.武汉理工大学高性能船舶技术教育部重点实验室，湖北武汉 430063； 2.武汉理工大学船海与能源动力工程学院，湖北武汉 430063； 3.国家知识产权局专利局专利审查协作湖北中心，湖北武汉 430063)

传统的船舶设计大多采用母型船改造法，该方法经常得到的是可行方案，而非最优方案。船型优化将计算流体力学与最优化理论结合并运用在船型设计过程中，可获得在一定约束条件下的性能最优的船型。

船体曲面修改是船型优化的基础。众多学者对船体曲面变形开展了深入研究，目前常用变形方法有自由变形方法、融合方法、径向基函数(RBF)变形方法以及参数化方法。李冬琴等[1]将非均匀有理B样条(NURBS)基函数和非均匀控制点框架结合，建立适用于船体几何的基于NURBS基函数的自由变形技术(NFFD)，以某CNG运输船为例完成球艏、船首和船尾部分的自动变形。冯佰威[2]在船体曲面NURBS表达的基础上，使用融合方法对DTMB5415声呐罩区域和1300TEU集装箱船球鼻艏区域进行了兴波阻力优化。Tahara等[3]使用融合方法将两种尾部型线差别较大的船型融合成新船型，避免了不符合实际的船型的出现。Kim等[4]将修改横剖面面积曲线方法和径向基函数插值方法相结合，达到对KCS船的船体曲面进行全局和局部修改的目的；沈通等[5]采用NURBS表达与RBF相结合的方法，实现了船体曲面的三维自动变形，研究结果表明，基于RBF的船体曲面修改方法适用于船型优化，符合工程实际需要；周昊等[6]采用RBF方法完成了DTMB5415船型的艏部型线优化。冯佰威等[7]采用RBF方法，完成了不同弗劳德数下的船舶阻力性能优化。张萍[8]提出了常规船型的参数化设计方法，由特征参数设计特征曲线，生成横剖面曲线，得到光顺的船体曲面。董枫[9]建立有效的船体参数化模型，针对复杂船型的各类横剖线构造了完整的形状参数集，采用改进的曲面放样法来构造船体曲面模型。胡春平[10]用Friendship软件对KCS船型进行全参数化建模，应用Shipflow软件以兴波阻力为目标对型线进行了优化。

文献[5]中对各种船体曲面变形方法进行了详细的对比分析。在诸多的船体曲面变形方法中，由于RBF船体曲面变形方法具有灵活高效，变形参数可控的特点，因此具有较强的工程应用价值。

在基于RBF的船体曲面变形方法中，支撑半径为关键参数，影响了船体曲面变形范围的大小。杨烜等[11]研究了径向基函数紧支撑集选取的问题，对任意标志点集，采用构造Delaunay三角剖分方法来确定基函数支撑集的大小；Zhang等[12]使用具有多个支撑半径的紧支撑径向基函数应用于心脑图像的配准中。

本团队在前期研究中采用Delaunay三角剖分方法求取支撑半径，并成功应用于船体曲面的光顺变形。但尚存在的问题是：在某些特殊情况下，Delaunay三角剖分无法获得支撑半径，限制了该方法的进一步推广应用。为解决上述问题，文中在前期研究基础上加以改进，提出了支撑半径的动态求取方法，并将该方法应用到标准船型优化中，以验证其可行性。

1 基于径向基函数插值的曲面变形方法

1.1 径向基函数及其插值技术

径向基函数是一种沿径向对称的标量函数，以空间中任一点Xi到某一中心x之间的欧氏距离r=‖X-Xi‖为自变量。各基函数的具体形式为[13]：

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

(1)

在船体曲面变形中应用RBF方法时，使用完整形式的插值方程：

(2)

船体曲面上中心点X=(x,y,z)移动的距离用s(X)表示，点的数量用n表示，Xi为任一点，p(X)为低阶多项式。其形式为

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

(3)

得到完整插值方程的关键在于求得λ和c。该方法的详细介绍请参考文献[13]。

1.2 支撑半径求取存在的问题

在船体曲面变形中，选取WendlandΨ3,1作为径向基函数插值的基函数，基函数选择参考文献[13]中的情况，该基函数形式如下：

(4)

在该函数中添加支撑半径R，使基函数变为如下形式：

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

(5)

支撑半径对基函数的影响范围进行控制，决定了曲面变形的质量。本团队主要采用Delaunay三角剖分的方法实现支撑半径的求取[14]。该方法尚存在的问题是：只考虑了可以进行三角剖分的情况，对于无法进行三角剖分的情况未作说明，且运用三角剖分求解支撑半径仅考虑到点的分布对支撑半径的影响，点的坐标变化量对支撑半径的影响并未考虑。这些缺点限制了船体曲面变形的应用范围。

2 支撑半径的动态求取

2.1 支撑半径求取分类

支撑半径动态求取的思路是：以可变点的数量为依据，要使变形后的船体曲面不出现鞍点且保持拓扑，需要综合考虑可变点的分布和可变点的坐标变化量对支撑半径的影响。在采用Delaunay三角剖分求取支撑半径时，重点在于可变点是否能构成三角形，着重考虑的是可变点的分布对支撑半径的影响。可变点的分布情况可分为图1中所示的几类，对于n<3或者n≥3且所有点共线，此时仅考虑可变点的分布情况，考虑无法采用Delaunay三角剖分时支撑半径的求取。同时将可变点坐标变化量考虑在内。

图1 可变点分布情况

2.1.1 考虑可变点坐标变化量的支撑半径求取

在船体曲面变形中要求保持拓扑。必要条件是函数S是连续的，并且雅可比矩阵的行列式必须在图像的每个点处为正：

det(∇s)>0

(6)

式中，∇s为雅可比矩阵。

为了证明变换的全局不变性属性，使用如下的定理[15]。

定理1 如果s:Ω→Rn，其中Ω是Rn中一个封闭的矩形区域，是可微分映射，使得雅可比矩阵∇s是Ω中所有x的P矩阵，那么s在Ω中是唯一的。

注意，如果n×n实矩阵的所有主子式都是正的，则称其为P矩阵。

为了变换函数s，分析了单个点p的雅可比矩阵，它包括一个恒等映射和一个局部基函数的叠加。对于2维的单个点有：

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

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

(7)

式中，φ是局部RBF基函数，Δ1和Δ2分别是在x和y坐标方向上从初始点p到目标点q的位移，并且p是任意位置的点。将式(7)代入式(6)，所得如下：

式中，φ表示φ(‖X-P‖)，r=‖X-P‖，如果令Δ=max(Δ1,Δ2)，在γ=π/4时，得到二维下的最坏情况：

(8)

对于三维的情况，获得了类似的条件：

如果再次设置Δ=max(Δ1,Δ2,Δ3)，获得在三维中最差的情况：

(9)

利用式(8)和式(9)，可以显示雅可比行列式的所有主子式分别对于二维和三维情况是正的。由此得出，如果式(8)和式(9)成立，则变换保留拓扑。∂φ/∂r的最小值取决于支撑半径R。表1总结了不同局部基函数的∂φ/∂r。在表1和式(8)、式(9)的前提下，可以导出R的条件。如表2所示。

表1 ∂φ/∂r的最小值

表2 给定位于Δ时最小的R

这些条件仅适用于在支撑半径R内没有其他点的孤立点。对于具有相交支撑区域的点，R的最小值还取决于点p的位置，这导致更复杂的计算。表2为选择支撑半径R的求取提供了线索。

在单标志点情况下，为了保持变换后的拓扑关系，WendlandΨ3,1基函数的支撑半径在三维时应满足R>3.66Δ，Δ=max(Δ1,Δ2,Δ3)，即支撑半径与点的坐标变化量有关，说明支撑半径有下界。

2.1.2 无法进行三角剖分时的支撑半径求取

无法进行三角剖分的情况如下：

(1)可变点的数量小于3。Delaunay三角剖分的基本思想是以三角网格的形式表征某一点与周围邻近点的拓扑关系，当可变点数量小于3个时无法构成三角形。

(2)可变点数量不小于3，但这些点共线。基础三角网增长法在寻找点和扩展边的时候要满足相应的准则，当不满足准则时无法进行三角剖分。

针对上述情况，给出的解决方法如下：

(1)可变点的数量为1时，仅考虑可变点坐标变化量的影响，选取的支撑半径为R=3.66Δ。

(2)可变点的数量不小于2时，此时为所有点共线的情况。对于WendlandΨ3,1基函数，杨烜等[11]已经证明：在二维两个控制点的情况下，只要保证两控制点间的区域不出现鞍点，就可以保证在其它位置不出现鞍点。并且当支撑半径取两点距离的2倍时，两点连线中点处的二阶导数小于0，不出现鞍点。记相邻两点间最大的欧氏距离为rd，仅考虑可变点分布的影响，取支撑半径为R=2rd。综合考虑可变点的分布和坐标变化量的影响，最终取支撑半径为R=max{2rd,3.66Δ}。

2.2 支撑半径的动态求取框架

支撑半径动态求取的具体算法框架如图2所示，详细的步骤可总结如下：

(1)首先在特定文件中读取出可变点的个数n(n为正整数)，然后求出各个点在每个方向(x,y,z)上的坐标变化量，并且求取出各变化量的最大值Δ。

(2)判断一用于确定可变点的数量n是否小于3。若n小于3，则因为点的个数太少无法构成三角形，从而不能进行Delaunay三角剖分，并且需要进行判断二；若n不小于3，则需要做进一步的判断，即判断三。

(3)判断二用于将可变点个数小于3的情况分为两种：当n等于1时，不用考虑可变点的分布对支撑半径的影响，支撑半径取为支撑半径的下界，即R=3.66Δ；当n等于2时，考虑可变点分布的影响，从而求取相邻两点间的最大距离rd，之后需要进行判断四。

(4)判断三用于将可变点个数不小于3的情况分为两种：当所有可变点共线的时候，不满足三角形的构成条件，无法进行Delaunay三角剖分，考虑可变点分布的影响，求取相邻可变点间的最大距离rd，之后进行判断四；当不是所有可变点都共线的时候，考虑可变点分布的影响，采用Delaunay三角剖分求取三角形的最大边长rd，之后进行判断四。

(5)判断四用于对可变点的分布位置和可变点的坐标变化量进行综合考量，同时满足变形之后不能产生鞍点，且保持拓扑结构，最终选取2rd和3.66Δ中较大的为支撑半径R。

图2 支撑半径动态算法框架

2.3 支撑半径的动态求取实例

将支撑半径的动态求取用于具体船型的变化，首先考虑可变点数量为2时的情况，该情况与可变点三点共线时求取策略相同。在DTMB5415船型的艏部选择两个可变点Y1、Y2，两点在设计水线附近，且坐标均在y轴方向变化，具体坐标数据见表3，经计算其3.66Δ为0.039 2，2rd为0.915 6，最终支撑半径为0.915 6；选择甲板边线以及声呐球艏处的点为固定点，如图3所示。运用支撑半径动态求取算法求取相应的支撑半径R，并且采用径向基函数插值方法得到变形后的船型，初始船和变形船的型线对比如图4所示。可以看出变形船在0水线(吃水为0处的水线)以下的型线未发生变化；在0水线以上的部分由横剖线图可知，其略微变肥，且首部型线内凹；半宽水线(首部)显示甲板边线没有发生变化，部分水线变凹。

表3 船型算例坐标数据表

图3 可变点数量为2时选点位置图

图4 可变点数量为2时型线对比图

考虑可变点数量为1时，在S60船型艏部选取一个可变点，可变点沿Y轴方向变化。可变点和固定点位置如图5所示，具体坐标数据如表4所示。运用支撑半径动态求取算法求取相应的支撑半径，由于可变点数量为1，因此支撑半径值为3.66Δ，经计算3.66Δ为0.023，因此其支撑半径为0.023。并且采用径向基函数插值方法得到变形后的船型，初始船和变形船的型线对比如图6所示，首部横剖线形状外凸，部分纵剖线呈S状。

图5 可变点数量为1时选点位置图

表4 船型算例坐标数据表

图6 可变点数量为1时型线对比图

综上所述，支撑半径动态求取方法可以进行特殊情况下支撑半径的求取，且所求取的支撑半径可以应用于径向基函数插值，实现船体曲面变形。

3 应用实例

3.1 优化问题描述

将改进的径向基插值曲面变形方法应用于SHIPMDO-WUT平台，对S60船型进行曲面的变形优化。用于变形优化的S60船模主尺度信息如表5所示。

表5 船模主尺度信息

优化目标定义如下(Fn为傅汝德数):

fobj=minCw,Fn=0.30

(10)

其中兴波阻力系数计算公式为

(11)

式中，ρ为海水的密度，V为航速，S为湿表面积，Rw为兴波阻力。

1)变量选择

选取如图7所示的4个曲面上的点，且每个点在y、z两个方向变化，例如点1的y坐标为设计变量Y1，点1的z坐标为设计变量Z1，依次类推，共4个可变点，8个坐标变化，即8个设计变量。具体设计变量及变量范围如表6所示。

图7 选点位置图

表6 设计变量及其范围

2)约束条件

静水力约束：优化船排水体积不得小于母型船。

型线约束：各主尺度要素船长、船宽、吃水均保持不变；首部轮廓线，基线保持不变，船体部分站位保持不变。固定点位置如图7所示。

3)优化算法

选择多目标粒子群算法(MOPSO)修改优化变量，设置20代，每代粒子数为50。在粒子群优化算法中，粒子速度和粒子位置更新见下式[16]。

(12)

式中:vi(t+1)为第t+1代、第i个粒子的速度;xi(t+1)为第t+1代、第i个粒子的位置；w为权重因子；c1、c2为学习因子；r1、r2表示0到1的随机数;α表示局部最优解;β表示全局最优解。当不满足优化收敛准则时，根据上式更新优化变量。

4)参数化变形模块

本流程的参数化变形模块即为基于径向基函数插值的船体曲面变形程序。该模块在迭代过程中读取更新的优化变量，通过径向基函数的插值操作，可实现船体曲面的光顺变形。

具体的优化流程如图8所示。

图8 优化流程图

3.2 优化结果分析

优化所得的1050组结果如图9所示。图中标注了初始船的兴波阻力系数位置，选择其中兴波阻力系数最小的船型为优化船。表7、表8分别列出了初始船和优化船的静水力数据、水动力数据。为了节省时间，优化过程中采用Shipflow软件计算船体兴波阻力系数；为了验证优化效果，优化结束后采用Fine/Marine软件分别计算初始船和优化船的总阻力。优化船相较初始船而言，排水体积增加0.83%，浮心纵向位置减小0.38%，湿表面积增加0.62%，兴波阻力减小10.46%，总阻力减小3.24%。

图9 优化收敛图

表7 静水力数据

表8 水动力数据

由式(11)可知，在排水体积和浮心纵向位置均满足约束的条件下，且航速一定时，优化船的湿表面积增加，兴波阻力系数减小，说明兴波阻力减小，且其减小是由船体型线的变化导致的；优化船总阻力减小，说明在优化时达到了通过兴波阻力系数减小从而使总阻力减小的目的。

图10显示的是初始船和优化船的型线对比图，图11为初始船和优化船的三维模型图。波形图和波切图的对比分别在图12、图13中显示。图13中LY为Y方向波切位置距船体中缝剖面的距离，图14为压力分布图。从型线对比图中明显可以看到型线变化主要在首部，产生了隐形球艏，部分纵剖面呈S型，对于高速时产生有利的兴波干扰有比较大的作用。由波切图可以看出，波高在波峰和波谷位置处的幅值减少，因此船首到船中部分(x/L为0～0.5，L为船长)的兴波幅值略微减小，船中到船尾部分(x/L为0.5～1.0)的兴波幅值明显减小，船尾及尾后0.5L部分(x/L为1.0～1.5)兴波幅值略微减小。说明隐形球艏的产生对船体兴波产生了有利干扰，使得兴波阻力减小。从压力分布图可以看出，优化船船首部分较初始船压力增加，船尾部分变化较小，因此其首尾压差阻力增加。这是由于船首产生了隐形球艏，湿表面积增加，黏性阻力增大。由于S60船型为中高速船，兴波阻力在总阻力中占比较大，而黏性阻力占比较小。因此优化船兴波阻力有明显降低，总阻力降低。