技术领域
[0001] 本
发明属于船舶性能计算领域，特别是一种船舶浮态计算方法。
背景技术
[0002] 船舶浮态计算是指计算船舶的船中
吃水TM、横倾
角θ及纵倾角 也可采用船舶艏吃水TF、船舶艉吃水TA及横倾角θ的表示方法。船舶任意装载状况下的浮态计算是船舶安全航行的重要前提，也是船舶稳性计算、船舶剪
力计算、船舶弯矩计算、船舶自动装载、船舶装载过程优化及确定船舶装卸货顺序的
基础。另外，船舶在装货过程中需要合理配载使船舶处于最佳纵倾状态下营运，以达到船舶节能的目的。
[0003] 目前船舶浮态计算的方法大体可以分为三类：
[0004] 一、第1类是依靠船舶静水力资料计算的“常规法”
[0005] 详见参考文献[1-5]，船舶大副进行手动计算时经常采用此种方法，设LBP为船舶两柱间长，XB为船舶
浮心纵向坐标，XF为水线面
漂心纵向坐标，MTC为每厘米纵倾力矩，为纵倾角，T为平均吃水，W为船舶总载重量，Pi和XGi分别为每个舱室的装货量及
重心纵向坐标，Pi由用户输入，XGi根据Pi查取舱容表插值获得，其基本步骤如下：
[0006] 1、由总载重量W查静水力表得到平均吃水T、XB、XF、MTC；
[0007] 2、计算重心纵向坐标： 式中n为所有重量
块个数。
[0008] 3、计算纵倾吃水差：
[0009] 4、计算纵倾角正切值：
[0010] 5、计算艏吃水：
[0011] 6、计算艉吃水：
[0012] 二、第2类是求解船舶空间浮态方程组的“矩阵法”
[0013] 船舶任意载况下的浮态平衡方程组为：
[0014]
[0015] 式中，Myz＝VXB-FD/ρ·XG为对中站面的体积矩，Mxy＝VZB-FD/ρ·ZG为对基平面的体积矩，Mxz＝VYB-FD/ρ·YG为对中线面体积矩，FD为船舶载重量。
[0016] 这是一个隐式非线性方程组，文献[6-10]提出采用逐次线性优化方法，通过船舶静水力原理，计算雅克比矩阵系数，将非线性方程组转换成线性方程组，这种方法称为“矩阵法”。
[0017] 三、第3类计算方法可归结为“优化法”
[0018] 主要有文献[11]提出的非线性规划法及文献[14-15]提出的遗传
算法。文献[11]把
优化算法中的非线性规划法用于浮态计算，采用惩罚函数法(SUMT)建立了船舶浮态约束优化模型：
[0019]
[0020] 构造惩罚函数如下：
[0021]
[0022] 式中f(x)为目标函数；r(k)、s(k)是在优化过程中随k增加而不断调整的变值参数，成为惩罚因子，k为
迭代次数。在优化计算过程中，取一系列r(k)、s(k)对惩罚函数进行一系列无约束最优化计算。当k→∞，惩罚函数的最优点逐步逼近原目标函数的最优点。
[0023] 文献[14-15]根据船舶自由漂浮的平衡条件，将自由浮态计算归结为多目标约束优化问题，并引入
遗传算法对该问题进行求解。
[0024] 目前这三类算法存在如下不足：
[0025] (1)第1类“常规法”存在如下假设：
[0026] 1)船舶在小角度纵倾时，其纵倾轴过初始水线面漂心的横轴，即漂心
位置不变；
[0027] 2)纵倾角较小时，
[0028] 3)纵稳性高度 和纵稳心半径R近似相等。
[0029] 从这些假设可以看出，常规法适合船舶小纵倾计算，纵倾较大时误差会增大，且不适用于船舶横倾较大的情况。
[0030] (2)第2类“矩阵法”在每一次迭代计算都要计算包含水线面面积、漂心、惯性矩、排
水体积及浮心等多项要素的雅克比矩阵，公式表达很繁琐，计算工作量大，若水线面为三体船或五体船等特种船舶时，程序处理困难。
[0031] (3)第3类“优化法”只需要计算船舶的排水体积和浮心，避免了雅克比系数的计算，减少了计算量。但是它们迭代次数较多，收敛速度缓慢，一般需要几百到几千次迭代才能找到优化结果，不适用于实时计算。
[0032] 本发明涉及的参考文献如下：
[0033] [1]盛振邦,刘应中.船舶原理(上册)[M].上海：上海交通大学出版社,2003:20-44。
[0034] [2]沈华,刘培学.船舶浮态的计算[J].大连海事大学学报,2004(03):18-20。
[0035] [3]王维宇,林洪波,何惠明.纵倾船舶浮态计算的新方法——三参数迭代法[J].交通部上海船舶运输科学研究所学报,1984(02):11-22。
[0036] [4]段兴锋,任鸿翔,东昉.纵倾状态船舶浮态的研究与计算[J].中国航海,2015,38(1)。
[0037] [5]刘春雷,尹勇,孙霄峰,等.船舶浮态计算的一种修正方法[J].大连海事大学学报,2014(04):1-6。
[0038] [6]赵晓非,王世连,李保拴.大纵倾下船舶浮态及初稳性计算[J].大连工学院学报,1982,01:77-85。
[0039] [7]赵晓非,王世连,李宝栓.船舶浮态计算方法及船舶浮态图谱[J].中国造船,1985(01):71-79。
[0040] [8]赵晓非,王世连,李保拴.大纵倾下船舶浮态及初稳性计算[J].大连工学院学报,1982(01):77-85。
[0041] [9]赵晓非,林焰.关于解船舶浮态问题的矩阵方法[J].中国造船,1985(03):57-66。
[0042] [10]林焰,李
铁骊,纪卓尚.破损船舶自由浮态计算[J].大连理工大学学报,2001(01):85-89。
[0043] [11]
马坤,张明霞,纪卓尚.基于非线性规划法的船舶浮态计算[J].大连理工大学学报,2003(03):329-331。
[0044] [12]孙承猛,刘寅东.船舶浮态计算的一种优化方法[J].大连海事大学学报,2006(02):39-41。
[0045] [13]胡丽芬,马坤,张凤香.破损舰船浮态稳性实时计算研究[J].武汉理工大学学报(交通科学与工程版),2009(06):1160-1163。
[0046] [14]陆丛红,林焰,纪卓尚.遗传算法在船舶自由浮态计算中的应用[J].上海交通大学学报,2005(05):701-705。
[0047] [15]金宁,谢田华,田恒斗,等.改进遗传算法在船舶自由浮态计算中的应用[J].中国航海,2007(01):10-12。
发明内容
[0048] 为解决
现有技术存在的上述问题，本发明要设计一种既可减少程序迭代次数，又能避免水线面要素求解的高
精度的基于STL模型的船舶任意浮态计算方法。
[0049] 为了实现上述目的，本发明的技术方案如下：一种基于STL模型的船舶任意浮态计算方法，包括如下步骤：
[0051] 采用两个左手坐标系：与静止海面相连接的固定坐标系ox'y'z'和与船体相连的船体坐标系oxyz，并且让两坐标系的原点重合，原点选在船舶基平面、中横剖面和中纵剖面的交点上。规定x轴指向船艏为正，y轴指向右舷为正，z轴向上为正。浮态参数选择符拉索夫参数，包括船中吃水TM、横倾角θ和纵倾角 并且规定向右舷倾斜的横倾角θ和向船艏倾斜的纵倾角 为正。
[0052] 两个坐标系之间的变换矩阵，从船体坐标系oxyz到固定坐标系ox'y'z'的坐标变换矩阵为：
[0053]
[0054] 从固定坐标系ox'y'z'到船体坐标系oxyz的坐标变换矩阵为：
[0055]
[0056] 式中，α为倾斜角，即倾斜水线面与基准面之间的夹角，
[0058] 船舶设计完成后，由船舶设计
软件导出船舶外壳及所有舱室的初始化图形交换规范模型即IGES模型，然后再通过3D建模软件转换成光
固化立体造型术模型即STL模型。
[0059] C、计算船舶装载量W及重心坐标(XG,YG,ZG)
[0060] 计算船舶总载重量及重心，总载重量包括船员及备品、货舱、压载舱、
淡水舱、燃油舱和
润滑油舱，总载重量及重心坐标按下式计算：
[0061]
[0062] 式中，W为船舶总载重量；W0为空船重量，XG0、YG0及ZG0分别为空船重心的纵向坐标、横向坐标及垂向坐标；Pi为各项装载重量；XGi、YGi及ZGi分别为各项装载重量重心的纵向坐标、横向坐标及垂向坐标。
[0063] D、按照“常规法”计算船舶船中吃水TM、横倾角正切值tanθ及纵倾角正切值计算流程如下：
[0064] D1、由总载重量W查静水力表得到平均吃水T、浮心纵向坐标XB、水线面漂心纵向坐标XF、每厘米纵倾力矩MTC和横稳心高度KMT，船舶静水力表由船舶装载手册获得。
[0065] D2、计算纵倾吃水差：
[0066] D3、计算纵倾角正切值：
[0067] D4、计算艏吃水：
[0068] D5、计算艉吃水：
[0069] D6、计算船中吃水：TM＝0.5(TA+TF)
[0070] D7、计算横倾角正切值：tanθ＝YG/(KMT-ZG)
[0071] 式中，Δ为船舶排水量，LBP为船舶两柱间长，XF为水线面漂心纵向坐标，XB为浮心纵向坐标，MTC为每厘米纵倾力矩，KMT为横稳心高度。
[0072] E、根据船中吃水、横倾角及纵倾角确定倾斜水线面的方程：
[0073]
[0074] F、通过倾斜水线面和船舶外壳STL模型求交，计算水线面下的船体排水体积V及浮心坐标(XB,YB,ZB)。
[0075] G、计算固定坐标系下的船舶重心坐标(X′G,Y′G,Z′G)及浮心坐标(X′B,Y′B,Z′B)[0076] 已知船体坐标系下重心坐标(XG,YG,ZG)、浮心坐标(XB,YB,ZB)计算固定坐标系下重心坐标(X′G,Y′G,Z′G)及浮心坐标(X′B,Y′B,Z′B)，其变换矩阵为：
[0077]
[0078]
[0080]
[0081]
[0082] H、判断浮态方程组是否满足平衡条件，即判断总载重量和排水量是否相等、重心和浮心是否在同一条垂线上。
[0083] H1、判断总载重量和排水量是否相等，即是否满足以下公式：
[0084] |W-V·ρ|＜ε1
[0085] 式中ε1为设定的迭代精度；
[0086] 如果不满足，则转步骤I；
[0087] H2、判断在固定坐标系下，重心与浮心的横向坐标是否相等，即是否满足以下公式：
[0088] |Y′G-Y′B|＜ε2
[0089] 式中ε2为设定的迭代精度；
[0090] 如果不满足，则转步骤I；
[0091] H3、在固定坐标系下，重心与浮心的纵向坐标是否相等，即是否满足以下公式：
[0092] |X′G-X′B|＜ε3
[0093] 式中ε3为设定的迭代精度；
[0094] 如果不满足，则转步骤I；
[0095] H4、输出当前浮态结果。
[0096] I、艏艉吃水及横倾角三参数迭代计算；
[0097] I1、进行第1次迭代；
[0098] I11、艏艉吃水及横倾角按照下式进行迭代计算：
[0099]
[0100] 式中，Aw为水线面面积，KML为纵稳心高度。
[0101] I12、计算当前浮态下的排水体积V1及浮心坐标(XB1,YB1,ZB1)
[0102] I13、按照步骤G的公式，计算固定坐标系下的浮心坐标(X′B1,Y′B1,Z′B1)[0103] I14、判断是否满足浮态方程组平衡条件精度要求：
[0104]
[0105] 若满足，输出当前浮态结果，程序结束；若不满足，令i＝2，转步骤I2继续迭代计算；
[0106] I2、进行第i次迭代；
[0107] I21、则第i次的迭代方程为
[0108]
[0109] I22、计算当前浮态下的排水体积V1及浮心坐标(XBi,YBi,ZBi)
[0110] I23、按照步骤G的公式，计算固定坐标系下的浮心坐标(X′Bi,Y′Bi,Z′Bi)[0111] I24、判断是否满足浮态方程组平衡条件精度要求：
[0112]
[0113] 若满足，输出当前浮态结果，程序结束；若不满足，令i＝i+1，转步骤I21继续迭代计算。
[0114] 与现有技术相比，本发明具有以下有益效果：
[0115] 1、本发明基于STL模型计算船舶任意装载状况下的浮态，STL(STereo Lithography)文件是一种3D模型文件格式，是3D SYSTEM公司于1988年制定的一个
接口协议。STL模型三维数据格式简单、跨平台性能优良、在计算几何特性时简单和与平面、曲面求交时计算速度快，广泛应用于科学计算
可视化、计算机动画、3D打印、
虚拟现实及快速成形等领域。
[0116] 2、相比于“常规法”只能计算船舶纵倾情况，且纵倾较大时误差会增大。本发明提出的三参数迭代方法，能同时计算船舶的横倾和纵倾，适用于船舶任意浮态，提高了计算精度。
[0117] 3、本发明相比于“矩阵法”减少了水线面面积、漂心、惯性矩等雅克比矩阵系数的求解，减少了计算量，只需计算出船舶任意浮态下的排水体积及浮心，程序健壮稳定。
[0118] 4、相比于“优化法”，本发明提出的方法收敛速度更快，对于纵倾较小的情况一般只需要1次迭代，对于大纵倾的情况一般需要3-5次迭代，提高了程序的实时性，具有一定的工程实用价值。
附图说明
[0119] 图1是本发明计算船舶浮态所包括的主要功能模块示意图。
[0121] 图3是本发明计算船舶浮态所采用的坐标系：船体坐标系和固定坐标系的横向视图。
[0122] 图4是本发明计算船舶浮态所采用的坐标系：船体坐标系和固定坐标系的纵向视图。
[0123] 图5是本发明所采用的散货船“SPRING COSMOS”外壳及所有舱室的STL模型。
[0124] 图6是本发明所采用的矿砂船“SHANDONG REN HE”外壳及所有舱室的STL模型。
[0125] 图7是本发明重心计算子程序流程图。
[0126] 图8是“常规法”计算船舶浮态流程图。
[0127] 图9是本发明静水力插值子程序流程图。
[0128] 图10是重心、浮心坐标系转换流程图。
[0129] 图11是15万次随机装载实验吃水差计算结果。
[0130] 图12是各吃水差下的平均迭代次数。
具体实施方式
[0131] 下面结合附图对本发明做进一步详细地描述，以散货船“SPRING COSMOS”及矿砂船“SHANDONG REN HE”为例：
[0132] 图1为本发明计算程序的主要功能模块，图2为船舶浮态计算的整体流程图，图3、4所示为程序采用的两个直角坐标系，本发明主要包括如下几个步骤：
[0133] 1、建立船舶外壳的STL(Stereo Lithography)模型。可由以下两个途径获得：(1)对于存在设计数据的船舶，可以船舶设计软件NAPA导出IGES(The Initial Graphics Exchange Specification)模型，然后再通过3D建模软件3Dmax转换成STL格式。(2)对于得不到船舶设计数据的船舶，就需要根据装载手册中的型值表、静水力数值表、总布置图、舱容图及舱容表等数据利用建模工具对外壳及舱室进行
三维重建。图5所示为散货船“SPRING COSMOS”外壳及所有舱室的STL模型，图6所示为矿砂船“SHANDONG REN HE”外壳的STL模型，这两船模型来源于船舶设计软件NAPA。
[0134] 2、计算船舶载重量及重心坐标，如图7所示。根据装载手册输入船舶空船重量及重心坐标，以船舶实际载况输入船上各个舱室的装载量及重心坐标。
[0135]
[0136] 按照上式计算出最终的总载重量W及重心坐标(XG,YG,ZG)。
[0137] 3、按照“常规法”计算船舶船中吃水TM、横倾角正切值tanθ及纵倾角正切值流程图如图8所示。由总载重量W查静水力表得到吃水T、浮心纵向坐标XB、水线面漂心纵向坐标XF、每厘米纵倾力矩MTC、横稳心高度KMT，流程图如图9所示。
[0138] 计算纵倾吃水差
[0139]
[0140] 计算纵倾角正切值
[0141] 计算艏吃水
[0142] 计算艉吃水
[0143] 计算船中吃水TM＝0.5(TA+TF)
[0144] 计算横倾角正切值，tanθ＝YG/(KMT-ZG)。
[0145] 最终得到当前载况的船中吃水、横倾角及纵倾角。
[0146] 4、根据船中吃水TM、横倾角正切值tanθ及纵倾角正切值 确定倾斜水线面的方程：
[0147]
[0148] 5、倾斜水线面和船舶外壳STL模型求交，计算水线面下的船体排水体积V及浮心坐标(XB,YB,ZB)。
[0149] 6、计算固定坐标系下的船舶重心坐标(X′G,Y′G,Z′G)及浮心坐标(X′B,Y′B,Z′B)，如图10所示。
[0150]
[0151]
[0152] 7、判断当前浮态是否满足浮态方程组平衡条件
[0153]
[0154] 式中，ε1为5吨，ε2为0.001m，ε3为0.001m，若满足，输出当前浮态计算值，程序结束。若有一个条件不满足按下式进行迭代
[0155]
[0156] 重复步骤4、5、6，直到精度满足下式
[0157]
[0158] 退出迭代，输出当前浮态计算结果，程序结束。
[0159] 本发明的收敛性试验验证如下：
[0160] 以38500散货船“DOLCE VITA”为例，在初始装载状况为空载的情况下进行了15万次随机装载实验，以最大舱容为限制在5个货舱随机产生装货量，如图11所示，吃水差的分布范围从-14m到14m。
[0161] 表1 15万次装载实验迭代次数
[0162]
[0163] 表1和图12所示为各个吃水差下的程序迭代次数，最大迭代次数为9次，出现在10～-14m吃水差，全部15万次装载方案的平均迭代次数3.23次。
[0164] 从图12可知，随着吃水差绝对值的减小，迭代次数逐渐减小：船舶吃水差位于-12～-14m平均迭代次数为6.48次；吃水差位于12～14m之间时平均迭代次数为6次；吃水差位于6～8m平均迭代次数为4.69次；当船舶吃水差处于-2～2m之间时，平均迭代次数为2次。船舶正常装货时，吃水差绝对值一般不会超过6m，一般不超过5次迭代即可收敛。
[0165] 实验结论：本发明
稳定性好，收敛速度较快，具有一定的工程应用价值。
[0166] 本发明的准确性实例验证：
[0167] 本发明以6.4万散货船“SPRING COSMOS”和25万吨矿砂船“SHANDONG REN HE”为例，进行了实例计算。表2所示为散货船“SPRING COSMOS”五种典型载况的浮态计算结果：LOAD00为空船载况、LOAD04为压载途中、LOAD08为压载到港、LOAD12为满载出港、LOAD23为
1、3、5舱隔舱装载出港。表3所示为矿砂船“SHANDONG REN HE”五种典型载况的浮态计算结果：LOAD00为空船载况、LOAD01为压载出港、LOAD05为压载到港、LOAD12为满载到港、LOAD21为2、4、6、8舱隔舱装载出港。
[0168] 分析表中数据可知：(1)“SPRING COSMOS”吃水差误差基本都在1cm左右，最大误差为1.8cm。“SHANDONG REN HE”吃水差误差基本都在1cm以下，最大误差为0.3cm。可见采用本发明计算的船舶浮态计算结果和装载手册给定值误差较小，具有较高的计算精度。(2)船舶空船载况下，船舶吃水差较大，本发明需要进行3次迭代，其他载况基本需要1次迭代即可获得较精确的结果，可见本发明收敛速度较快，程序实时性好，具有一定的工程实用价值。
[0169] 表2 6.4W吨散货船“SPRING COSMOS”浮态计算结果
[0170]
[0171] 表3 25W吨矿砂船“SHANDONG REN HE”浮态计算结果
[0172]
