NURBS流曲线造型新方法及其在船舶设计中应用

张 彦 儒, 林 焰,2, 陆 丛 红*, 纪 卓 尚

( 1.大连理工大学 船舶工程学院, 辽宁 大连 116024； 2.大连理工大学 工业装备结构分析国家重点实验室, 辽宁 大连 116024 )

NURBS流曲线造型新方法及其在船舶设计中应用

张 彦 儒1, 林 焰1,2, 陆 丛 红*1, 纪 卓 尚1

( 1.大连理工大学 船舶工程学院, 辽宁 大连 116024； 2.大连理工大学 工业装备结构分析国家重点实验室, 辽宁 大连 116024 )

针对当前流曲线造型方法仅被应用于已知型值点及其切矢的问题，在其他约束(面积、面积心、流场特征值)问题中的应用，以及与NURBS设计方法的结合仍有待发展研究的现状，将NURBS与流曲线造型技术相结合，提出了NURBS流曲线单元的概念以及在其他几何约束下，应用NURBS流曲线造型的新方法．该方法将待设计流曲线在已知切矢的型值点处分段以确定NURBS流曲线单元的个数，将NURBS流曲线单元控制参数的并集作为设计变量，将几何约束处理为曲线设计的目标函数，建立普适性的NURBS流曲线造型框架．该框架模型可用进化算法进行求解．通过设计实例，验证了该方法在相关曲线设计方面的可行性和有效性．

船体曲面；流曲线；NURBS；功能曲线；几何约束

0 引 言

在运动物体(如汽车、飞机、船舶等)外形设计中，运动物体在空气、水流等流体中作相对运动，流体对其产生阻力，因此其外形设计变得十分重要．运动物体外形的光滑与否将直接影响其运动性能．而具有流线形外形的运动物体不仅外观漂亮，并且能极大地减少前进过程中流体对物体的阻力．针对运动物体的外形设计，以流体力学为背景的流曲线曲面的造型方法被提出，这种方法称为流曲线曲面造型技术．纵观现有文献，可以将流曲线曲面造型技术分为两类：第一类是用切矢的积分曲线进行曲线曲面造型的技术[1-6]，其本质是曲线切矢的拼接，是以弧长为参数的三维插值曲线．国内研究者在国外研究的基础上，在旋转曲面[7]、非均匀细分[8]、双三次非均匀B样条曲面的G1混合[9]等方面进行了研究．但是这一类流曲线的切矢并非流场中的速度矢量，因此本文定义该类方法为广义的流曲线曲面，它并非是流线形设计中的功能曲线．第二类是基于流体力学原理，将流体力学中流函数的概念引入CAD来构造流曲线曲面．1996年卢小林等给出了流曲线曲面的概念及多项式和B网构造流曲线曲面的基本方法[10]，1997年游世辉等给出了流体曲面设计的一般流体光顺条件[11]．2006年刘国伟提出了来流为直匀流的流曲线造型方法，并于2010年应用于高速列车头部曲线设计[12-13]．2017年徐阳等提出一种面向叶型和翼型的考虑流线场约束的NURBS曲线拟合方法[14]．上述文献的造型方法，仅局限在流曲线上型值点及其切矢的设计意图；而流曲线设计，作为一种功能曲线设计，不只是考虑曲线上型值点和切矢的几何条件，有时还需要考虑曲线所围面积、面积心等几何约束，但未见相关报道．另外，1991年国际标准化组织(ISO)颁布的工业产品数据交换标准STEP中，把NURBS作为定义工业产品几何形状的唯一数学方法．因此，研究如何将NURBS曲线与流曲线相结合，对流曲线造型的发展将会有非常积极的推动作用．并且，随着计算机技术、数值计算方法和优化方法的发展，船舶的外形设计由常规的正向设计逐渐走向反向设计[15]，因此研究如何对周围流场信息进行提取，反向指导功能曲线设计具有现实的工程意义．总体来说，流曲线曲面造型技术的基本理论、数学模型和相关算法仍然不够完善，还有待进一步研究扩展．

本文在第二类流曲线造型技术的基础上，将待设计流曲线用NURBS流曲线单元组合设计，并给出在曲线几何约束(型值点、面积、面积心、流场中特征点特征值)条件下，NURBS流曲线的数值造型方法，以弥补流曲线造型只能处理型值点、切矢等的造型方法的不足，并可应用于反向设计．

1 流曲线数值求解模型

第二类流曲线基本是以理想不可压缩流体平面定常无旋运动为力学模型构造柔性样条，将满足流函数的全部或部分流线作为待求的功能性曲线．它是具有流线特征、可描述物理上可能的流动、可携带流场物理信息的功能曲线[12]．

由流体力学可知二维无旋不可压缩定常流动的流函数ψ(x,y)满足拉普拉斯方程：

(1)

流函数的等值线ψ(x,y)=C(C为常数)就是流场中的流线，流线上点的切矢与流场中的速度矢量相一致．如图1所示，同一个物体在不同的流场a和b中运动有不同的流线，并且即使是在同一个流场中，常数C取不同的值也对应着弯曲程度不同的流线．因此基于拉普拉斯方程构造的流曲线具有良好的可塑性．另外，由物面的不可渗透边界条件及流线的流体不可穿越性可知，流场中的任意流线均可用同位置、同形状的物面代替．因此流曲线造型的实质是求解满足给定流动状态的拉普拉斯方程以得到在特定位置、满足一定形状约束的流线．

(a) 流场a

(b) 流场b

求解拉普拉斯方程，一种是给定物体的外形和绕流形式的前提下求解流函数，这一类称之为正问题，可以应用常规的求解偏微分方程的方法进行求解．而另一种，是在给定部分流函数的前提下，求解物体形状，这一类称之为反问题．由前面叙述可知，流曲线造型属于反问题．

1.1 流曲线造型的数值解法

流曲线造型数值解法的实质是拉普拉斯方程第一边值问题的逆问题求解．因此先给出其求解公式，然后给出几何约束条件的处理，以及需要注意的几个补充方程．本文采用五点差分格式，对拉普拉斯方程进行求解．五点差分格式为

ψi,j=(ψi,j+1+ψi-1,j+ψi+1,j+ψi,j-1)/4

(2)

1.1.1 计算域 如图2所示，给出来流为直匀流条件下，坐标系和边界条件的定义．曲线OEF为待求流曲线．AB、BC、CD、AO、FD为固定边界，OF为控制边界，流曲线OEF的形状完全由控制边界OF的边界值确定．h为x、y方向上的网格步长．○点称为网格内点，×点称为边界点．x轴定义在计算域边界AD上，指向固定边界CD为正．y轴指向固定边界BC为正．坐标原点由待求流曲线OEF的几何特点决定．设来流方向平行于x轴，沿x轴正向．

图2 计算域

文献[13]未给出固定边界的数值确定方法．经过反复比较和分析，证实如果在曲线设计阶段，给出了流场中网格点值或者流场中特定点切线向量时，将进流段AO长度和去流段FD长度，以及计算域宽度AB按给定的流场区域进行计算即可；否则可以分别按图2所示，相对OF长度的对应比例选取．

1.1.2 边界条件 设来流为直匀流，平行于x轴，沿x轴正向，速度大小为vs．因固定边界处流场信息与待求流曲线形状无关，所以固定边界处边界条件可按直匀流流场特征给出．边界条件可定义为

(3)

其中yB表示边界点B的y坐标．从图2和式(3)可知，边界条件中只有控制边界OF段为未知．文献[13]给出了已知零流曲线(物面)上部分型值点及其切矢的条件下，将反问题转化为正问题，运用五点差分法求解拉普拉斯方程的数值造型方法．该方法虽然给出了切矢转化为网格点上流函数值的方法，但并未给出如何在满足计算精度的网格步长下，使型值点落到网格点上，而且未知控制边界上的边界点的个数只能等于已知内点个数，因此其网格不能细分，计算精度不高．本文不考虑控制边界OF段的边界值，直接将待求零流曲线OEF用NURBS表示，进行流曲线造型设计．给出补充方程：

ψ(x,y)=0； (x,y)∈OEF

(4)

这样，只要零流曲线给定，拉普拉斯方程就可以求解．考虑到待求的零流曲线作为控制边界时，边界可能不通过网格点，在给定特征点后划分网格，边界点也可能不通过网格点，给出上述情况的临近边界点网格点的计算方法如下．如图3所示在网格点(i,j)处的导数，可按下述方法计算．

图3 边界不通过网格点

根据Taylor展开公式：

(5)

(6)

消去一阶导数∂ψi,j/∂x，可以得到网格点(i,j)处关于x的二阶导数：

(7)

同理，可以得到网格点(i,j)处关于y的二阶导数：

(8)

因此拉普拉斯方法在网格点(i,j)处可近似为

(9)

其中ψA和ψB为拉普拉斯方程的边界值，在本文中因为是零流曲线，所以取值为0．

1.2 零流曲线的NURBS表示

一条p次NURBS曲线定义为

(10)

经过推导，可以得到NURBS曲线在端点处(u=0，u=1)的一阶导矢，其表达式如下：

(11)

考虑到设计流曲线一般不会发生折转，经过反复比较和分析，证实了利用三次NURBS曲线及5个控制顶点来设计一般自由曲线，可以满足工程精度和灵活修改等要求．将这个自由曲线段称为NURBS流曲线单元．

如图4所示P0是流曲线单元的起点，P4是流曲线单元的终点．假设xi、yi和wi分别是控制顶点Pi(i=0，1，…，4)的纵向坐标、横向坐标及该控制顶点对应的权因子．自由段曲线的首尾控制点权因子w0和w4可设置为1．节点矢量U可以设置为(0 0 0 0 0.5 1 1 1 1)．首末控制顶点P0和P4为已知物面(零流曲线)的型值点，y1和y3可以由式(11)以及w1和w3联合确定，因此流曲线单元的设计变量为x1、x2、x3、y2、w1、w2、w3．

图4 流曲线单元控制顶点分布

将已知的零流曲线上的型值点分为两种：一种为型值点和型值点处曲线切矢已知，另一种为型值点已知而型值点处切矢未知．用第一种型值点将零流曲线进行分段，每一段用一个流曲线单元表示．第二种型值点作为几何约束进行处理．

1.3 几何约束的处理

本文将几何约束转化为设计目标进行处理，即将设计曲线相关几何量与几何约束间的相对误差的线性加权值最小作为设计目标．即

(12)

式中：n为几何约束的个数，αi为各个几何约束的权重系数，Fi为几何约束的目标函数．下面分别叙述两种典型的几何约束处理方式．

1.3.1 流场中特征值的处理 流曲线造型，不能只考虑物面(零流曲线)上的切矢处理，更应该考虑周边流场信息的提取，并依此来指导物面(零流曲线)的设计．流场中的特征值一般分为特征点的速度矢量和特征点的流函数两类．由流线的性质可知：流线上任意一点的速度矢量在该点必与流线相切．如图5所示，流场中网格点(i,j)处的切矢就是实际流场中该点的速度矢量v．因此可以将特征点上的速度矢量按如下方式转化为特征点流函数值进行处理．

vx、vy与流函数之间的关系及差分方程为

(13)

图5 特征点上速度矢量处理

网格点(i,j)处的流函数值ψi,j与网格步长h已知，当vx、vy已知时，可以由式(13)推出相邻网格点的流函数值．这样，可以将速度矢量几何约束的设计目标转化为最小化对应网格点处流函数值的相对误差的最大值，设ψdi,j+1和ψdi+1,j为通过本文方法设计的流曲线流场中网格点(i,j+1)和(i+1,j)的流函数值，则设计目标可以表示为

minF1=max{|ψdi,j+1-ψi,j+1|/|ψi,j+1|,

|ψdi+1,j-ψi+1,j|/|ψi+1,j|,

|ψdi,j-ψi,j|/|ψi,j|}

(14)

1.3.2 未知切矢型值点的处理 假设欲使设计曲线通过型值点列Pod=(Pod1Pod2…Podn)，流曲线上对应型值点x坐标的y坐标为(Pd1.yPd2.y…Pdn.y)．因此设计目标为最小化各型值点y坐标值与对应的流曲线上点的y坐标值的相对误差中的最大值，即

minF2=max{|Pd1.y-Pod1.y|/|Pod1.y|,

|Pd2.y-Pod2.y|/|Pod2.y|,…,

|Pdn.y-Podn.y|/|Podn.y|}

(15)

1.3.3 流曲线面积和面积心的处理 设Awo为所设计的曲线和控制边界围成的区域要求的面积，(Laox,Laoy)为该区域面积的面积心坐标，Awd和(Ladx,Lady)分别为通过本文设计方法得到的流曲线的相应参数，目标函数设为最小化设计曲线对应区域的面积相对误差和面积心坐标的相对误差之和，即

minF3=min{|Awd-Awo|/|Awo|+

|Ladx-Laox|/|Laox|+

|Lady-Laoy|/|Laoy|}

(16)

2 流曲线造型步骤

本文采用以下方法做NURBS流曲线造型：

Step1研究待设计流曲线的几何性质特征，给出最少控制顶点的NURBS设计模型即流曲线单元．

Step2给出已知切矢的型值点．然后用已知切矢的型值点将待设计流曲线进行分段，分别对每一段用流曲线单元进行设计，最后汇总给出总的设计变量．

Step3给出其他几何约束条件，根据1.3约束处理方法将几何约束转化为设计目标．

Step4选用进化算法，将Step2中确定的设计变量和Step3中确定的目标函数代入进化算法进行求解．

Step5输出设计得到的流曲线及其周边流场信息．

3 数值算例

根据本文的造型方法，采用人机交互的遗传算法[16]对7 000 t 散货船(设计航速18 kn)设计水线进行方案设计的数值算例．为了配合图2计算域的坐标，本船坐标系设置为坐标系原点在舯横剖面与中纵剖线的交点，x轴指向船尾为正，y轴指向左舷为正．按水线特征将水线分为进流段和去流段．表1列出了已知条件和设计目标．其中Awo表示水线进流段或去流段与船中纵剖线和所设计曲线自由段的端点处横剖面所围面积，(Laox,Laoy)表示面积心坐标，Pst表示流曲线单元起点坐标，vst表示起点处切矢，Ped表示流曲线单元终点坐标，ved表示终点处切矢，ψC为计算流场中特征点处的流函数值．

图6给出了本文方法得到的设计水线结果以及控制顶点分布，图7给出了周围流场信息．

表1 设计水线的已知条件和设计目标

Tab.1 The known conditions and design objective of the design waterline

已知条件Pst/mPed/mvstved进流段(-50．764，0．405)(-23．257，8．600)(79．4847．12)(397．5 0)去流段(25．5，8．6)(53．799，3．348)(322．8 0)(1013．4-262)设计目标Awo/m2(Laox，Laoy)/m特征点/mψC进流段168．128(-33．958，3．557)(-50．995，1)4．407去流段187．359(37．656，3．515)(53．055，4)10．255

图6 设计水线的结果

图7 周围流场信息

表2列出了本文方法得到的水线控制顶点坐标．其中，x和y分别是控制顶点的x、y坐标，w是对应控制点的权值．表3列出了算法得到的水线的特征参数以及与设计目标的相对误差|Er|．由于进流段端部存在拐点，相比去流段的误差较大．但对于各项几何设计目标，已基本达到了工程上的初始设计精度要求(工程精度要求一般为0.5%)．

表2 设计水线控制顶点坐标

表3 设计水线设计结果与设计目标间的差异

4 结 语

本文将流曲线造型方法进行扩展，结合NURBS和进化计算方法，给出了在给定约束条件下，利用拉普拉斯方程第一边值问题的逆问题数值解进行流曲线造型的新方法．研究结果表明：本文可以运用周围流场的信息来指导流曲线的设计，并且可以在给定来流的条件下初步对其进行评估．

本文对于流曲线造型方法的扩展仅做了初步探讨，在以下方面仍需进一步研究：如何将二维设计转化为三维设计；来流为非直匀流时，流曲线的造型方法；如何提取周围三维流场信息以指导二维特征线设计；采用并行计算方法以加快流曲线设计速度．总之，流曲线曲面造型技术有待进一步研究，应用领域仍有待扩展．

另外，算例中，进流段误差相比去流段误差较大，是因为设计目标水线的进流段端部存在拐点，逼近较困难，造成设计结果与设计目标间面积、面积心及流函数值的误差较大．在五点差分方法中，对计算域进行网格离散后，网格点值对边界值比较敏感，微小的边界值波动，也会引起较大的流函数值的误差．因此流曲线设计结果取决于相关几何特性的误差以及设计结果形状是否满足设计预期．而流函数值的作用是作为指导约束，与相关几何特性一起指导流曲线的设计，使其设计形状逼近设计预期．其允许误差与网格点所在位置和设计目标形状有关，它的具体规律将是后续的研究内容．

[1] MIURA K T, WANG L, CHENG F. Streamline modeling with subdivision surfaces on the Gaussian sphere [J].ComputerAidedDesign, 2001,33(13): 975-987.

[2] SUZUKI Y, TANAKA I, MIURA K T,etal. Streamline modeling based on potential flow [C] //Proceedings-9thPacificConferenceonComputerGraphicsandApplications,PacificGraphics2001. Tokyo: IEEE Computer Society, 2001:305-312.

[3] MIURA K T, CHENG Fuhua, WANG Lazhu. Fine tuning: curve and surface deformation by scaling derivatives [C] //Proceedings-9thPacificConferenceonComputerGraphicsandApplications,PacificGraphics2001. Tokyo: IEEE Computer Society, 2001: 150-159.

[4] MIURA K T, SAKIYAMA N, KANEKO T. Streamline modeling-Designing free-form surfaces by blending tangent vectors [J].IFIPAdvancesinInformationandCommunicationTechnology, 2001,75:235-248.

[5] SAKIYAMA N, MIURA K T, TAKAHASHI T,etal. Fair curve and surface design system using tangent control [C] //PG′98Proceedingsofthe6thPacificConferenceonComputerGraphicsandApplications. Washington D C: IEEE Computer Society, 1998: 27-36.

[6] MIURA K T. Unit quaternion integral curve: A new type of fair free-form curves [J].ComputerAidedGeometricDesign, 2000,17(1): 39-58.

[7] 李秀娟,刘 浩,何 钢,等. 基于流曲线的旋转曲面的生成和形状调整[J]. 机械科学与技术, 2008,27(3):326-329.

LI Xiujuan, LIU Hao, HE Gang,etal. Generation and shape adjustment of revolution surface based on stream curve [J].MechanicalScienceandTechnologyforAerospaceEngineering, 2008,27(3):326-329. (in Chinese)

[8] 李秀娟,刘 浩,何 钢,等. 基于非均匀细分的流曲线曲面[J]. 机械科学与技术, 2008,27(10):1216-1219, 1224.

LI Xiujuan, LIU Hao, HE Gang,etal. Stream curve and surface based on non-uniform subdivision [J].MechanicalScienceandTechnologyforAerospaceEngineering, 2008,27(10): 1216-1219, 1224. (in Chinese)

[9] 李秀娟,廖文和,刘 浩,等. 基于流曲线曲面的双三次非均匀B样条曲面G1混合[J]. 中国科学院大学学报, 2008,25(4):554-559.

LI Xiujuan, LIAO Wenhe, LIU Hao,etal.G1continuity blending between bicubic non-uniform B-spline surfaces with stream curve and surface [J].JournaloftheGraduateSchooloftheChineseAcademyofSciences, 2008,25(4):554-559. (in Chinese)

[10] 卢小林,马利庄,何志均. 流曲线曲面的概念和构造方法[J]. 计算机工程, 1996,22(3):26-31.

LU Xiaolin, MA Lizhuang, HE Zhijun. The conception and construction method for fluid curve [J].ComputerEngineering, 1996,22(3):26-31. (in Chinese)

[11] 游世辉,卿启湘,卿 钧. 计算机辅助流体曲面优化设计的探讨[J]. 湖南大学学报, 1997,24(1):71-75.

YOU Shihui, QING Qixiang, QING Jun. Research of computer aided optimal design for fluid surface [J].JournalofHunanUniversity, 1997,24(1):71-75. (in Chinese)

[12] 刘国伟. 流曲线概念及造型方法[J]. 计算机辅助设计与图形学学报, 2006,18(12):1891-1896.

LIU Guowei. Conception and modeling method of stream curve [J].JournalofComputer-AidedDesign&ComputerGraphics, 2006,18(12):1891-1896. (in Chinese)

[13] LIU Guowei, LIU Weixiao. Stream curve and numerical modeling for high-speed train [C] //InternationalConferenceonOptoelectronicsandImageProcessing. Piscataway:IEEE Computer Society, 2010: 407-410.

[14] 徐 阳,刘 强. 考虑流线场约束的NURBS曲线拟合方法[J]. 计算机辅助设计与图形学学报, 2017,29(1):137-144.

XU Yang, LIU Qiang. NURBS curve fitting method considering constraint of streamline field [J].JournalofComputer-AidedDesign&ComputerGraphics, 2017,29(1):137-144. (in Chinese)

[15] 赵 峰,李胜忠,杨 磊,等. 基于CFD的船型优化设计研究进展综述[J]. 船舶力学, 2010,14(7):812-820.

ZHAO Feng, LI Shengzhong, YANG Lei,etal. An overview on the design optimization of ship hull based on CFD techniques [J].JournalofShipMechanics, 2010,14(7):812-820. (in Chinese)

[16] 王运龙,王 晨,彭 飞,等. 基于人机结合遗传算法的船舶管路三维布局优化设计[J]. 中国造船, 2015,56(1):196-202.

WANG Yunlong, WANG Chen, PENG Fei,etal. Intelligent layout design of ship pipes based on genetic algorithm with human-computer cooperation [J].ShipbuildingofChina, 2015,56(1):196-202. (in Chinese)

NewNURBSstreamcurvemodelingmethodanditsapplicationtoshipdesign

ZHANGYanru1,LINYan1,2,LUConghong*1,JIZhuoshang1

(1.SchoolofNavalArchitecture&OceanEngineering,DalianUniversityofTechnology,Dalian116024,China；2.StateKeyLaboratoryofStructuralAnalysisforIndustrialEquipment,DalianUniversityofTechnology,Dalian116024,China)

The stream curve modeling method has only been applied to the problems with data points and their tangent vectors being given. The application to the problems with other constraints such as area, centroid of area, the eigenvalue of the flow field, etc. and the integration with NURBS need more research and development. NURBS and stream curve modeling method are integrated and the concept of NURBS stream curve unit is proposed. Furthermore, a new NURBS stream curve modeling method based on other geometric constraints is proposed. In this method, the stream curve is segmented at the data points with the given tangent vectors, and the number of the NURBS stream curve cell is determined. The union sets of the controls parameter of the NURBS stream curve cell are designated as the design variables. The geometric constraints are regarded as the design objective function. The universal framework of NURBS stream curve modeling is built. This framework model can be solved by an evolutionary algorithm. The design instances indicate that this method is effective and feasible in the related curve design.

hull surface; stream curve; NURBS; functional curve; geometric constraint

1000-8608(2017)06-0564-07

U662.2

A

10.7511/dllgxb201706003

2017-04-05；

2017-09-25．

国家自然科学基金资助项目(E091002-51109033).

张彦儒(1985-)，男，博士生，E-mail:03401011@163.com；林 焰(1963-)，男，博士，教授，博士生导师，E-mail:linyanly@dlut.edu.cn；陆丛红*(1972-)，女，博士，副教授，E-mail:lchcad@dlut.edu.cn.