一种液货船配载仪中船舶破损GZ曲线计算方法

一种液货船配载仪中船舶破损gz曲线计算方法
技术领域
1.本发明属于船舶性能计算领域，特别是一种液货船配载仪中船舶破损gz曲线计算方法。


背景技术：

2.如图1所示，船舶gz曲线是指船舶某装载状况下在各个横倾角度下(0
°
、5
°
、10
°
、20
°
、30
°
、40
°
、50
°
、60
°
及70
°
)船舶重力作用线与浮力作用线之间的距离所连成的曲线。船舶gz曲线是船舶破舱稳性计算的重要组成部分，是船舶安全航行的重要前提。
3.目前船舶gz曲线计算的方法大体包括以下类别：
4.第一类是依靠船舶静水力资料计算的“离线插值法”，船舶大副进行手动计算时经常采用此种方法，该方法的基本路线是：设xb为船舶浮心纵向坐标，xf为水线面漂心纵向坐标，mtc为每厘米纵倾力矩，为纵倾角，t为船中吃水，w为船舶总载重量，pi和x
gi
分别为每个舱室的装货量及重心纵向坐标，pi由用户输入，x
gi
根据pi查取舱容表插值获得，xg及zg为当前载况的总重心纵向及垂向坐标，其基本步骤如下：
5.1、由总载重量w查静水力表得到船中吃水t、xb、xf、mtc；
6.2、计算重心纵向坐标：式中n为所有重量块个数。
7.3、计算纵倾吃水差：式中δ为船舶排水量。
8.4、根据吃水差及船中吃水查询装载手册中稳性横截力臂表值kn
9.5、计算gz值：gz＝kn-zggsinθ
10.第二类是“实时计算法”，设船舶破损载况下的平衡方程组为：
[0011][0012]
式中，m
xy
＝vz
b-fd/ρ
·
zg为对基平面的体积矩，m
xz
＝vy
b-fd/ρ
·
yg为对中线面体积矩，fd为船舶载重量。
[0013]
这是一个隐式非线性方程组，采用逐次线性优化方法，通过船舶静水力原理，计算雅克比矩阵系数，将非线性方程组转换成线性方程组，这种方法称为“实时计算法”。
[0014]
第三类计算方法可归结为“优化法”，主要有非线性规划法及遗传算法。
[0015]
目前这三类算法存在如下不足：
[0016]
(1)第一类“离线插值法”存在如下不足：离线插值法依赖于船舶装载手册，装载手册中一般会给出尾倾5m到首倾3m的稳性力臂数据，在实际装载工况中当船舶吃水差不在-5m～3m的范围时，计算误差会增大。
[0017]
(2)第二类“实时法”在每一次迭代计算都要计算包含水线面面积、漂心、惯性矩、排水体积及浮心等多项要素的雅克比系数矩阵，公式表达很繁琐，计算工作量大，需要实时
计算船舶水线面。
[0018]
(3)第三类“优化法”只需要计算船舶的排水体积和浮心，避免了雅克比系数矩阵的计算，减少了计算量。但是它们迭代次数较多，收敛速度缓慢，一般需要几百到几千次迭代才能找到优化结果，实时性较差。
[0019]
本发明涉及的参考文献如下：
[0020]
[1]盛振邦，刘应中.船舶原理(上册)[m].上海：上海交通大学出版社，2003:20-44。
[0021]
[2]王智洲，孙霄峰，尹勇，等.散货船配载仪中完整稳性计算方法研究[j].大连理工大学学报，2017，57(5):6.
[0022]
[3]马坤，李作志，杨洋.舰船完整稳性实时计算方法研究[j].大连理工大学学报，2007，47(3):4.
[0023]
[4]崔阳.基于catia三维建模的船舶稳性计算方法研究[d].大连理工大学，2015.
[0024]
[5]令波，张正艺，解德.基于abaqus的破损船舶静稳性曲线直接计算法研究[j].中国水运，2020(12):2.
[0025]
[6]沈晔，孙霄峰，尹勇，等.散货船配载仪中破舱稳性计算[j].大连理工大学学报，2018，58(1):7.
[0026]
[7]肖维维，李俊华，何刚，等.船舶破舱稳性计算方法研究[j].中国水运(理论版)，2006(11):12-14.
[0027]
[8]蒋昌师，凌伟.小水线面双体船破舱稳性计算分析[c]//北京造船工程学会2016-2017年学术论文集..
[0028]
[9]黄武刚.基于foran的船舶破舱稳性计算分析研究[c]//全国船舶稳性学术研讨会.2014.
[0029]
[10]张明霞.基于nurbs曲面的船舶破舱稳性计算方法研究[d].大连理工大学，2002.
[0030]
[11]赵晓非，王世连，李保拴.大纵倾下船舶浮态及初稳性计算[j].大连工学院学报，1982，01:77-85。
[0031]
[12]林焰，李铁骊，纪卓尚.破损船舶自由浮态计算[j].大连理工大学学报，2001(01):85-89。
[0032]
[13]刘春雷，尹勇，孙霄峰，等.基于stl模型的散货船破损浮态计算[j].系统仿真学报，2016，28(10):8.
[0033]
[14]姜玲.船舶稳性三维计算方法研究及软件开发[d].大连海事大学，2015.
[0034]
[15]胡丽芬，马坤，张凤香.破损舰船浮态稳性实时计算研究[j].武汉理工大学学报(交通科学与工程版)，2009(06):1160-1163。
[0035]
[16]胡丽芬.舰船抗沉辅助决策系统研究.大连理工大学，2010.
[0036]
[17]陆丛红，林焰，纪卓尚.遗传算法在船舶自由浮态计算中的应用[j].上海交通大学学报，2005，39(5):6.


技术实现要素：

[0037]
为解决现有技术存在的上述问题，本发明要设计一种既可减少程序迭代次数，又
能避免水线面要素求解的液货船配载仪中船舶破损gz曲线计算方法。
[0038]
为了实现上述目的，本发明的技术方案如下：一种液货船配载仪中船舶破损gz曲线计算方法，包括如下步骤：
[0039]
a、建立坐标系
[0040]
采用左手坐标系，建立船体坐标系oxyz，原点选在船舶基平面、中横剖面和中纵剖面的交点上。规定x轴指向船艏为正，y轴指向右舷为正，z轴向上为正。
[0041]
b、建立船舶平衡方程组
[0042]
船舶破损进水后的稳定状态需要满足2个条件，即：重力和浮力相等、重心和浮心在一条直线上，具体公式如下：
[0043][0044]
式中，ρ为海水密度，v为船舶排水体积，δ为船舶排水量，xb、yb、zb分别为船舶浮心纵向、横向及垂向坐标，xg、yg、zg分别为船舶重心纵向、横向及垂向坐标，θ为船舶横倾角，为船舶纵倾角。
[0045]
c、简化雅克比系数矩阵
[0046]
对于步骤b建立的如下平衡方程组：
[0047][0048]
引入向量表示如下：
[0049][0050]
式中，t为船中吃水。
[0051]
使用牛顿迭代法得到线性化方程：
[0052][0053]
式中，δt为船中吃水改变量、为船舶纵倾角正切值改变量，为雅克比系数矩阵。
[0054]
根据船舶静力学原理，各偏导数计算公式为：
[0055][0056]
则最终雅可比系数矩阵为：
[0057][0058]
式中，s为船舶倾斜水线面在基平面的投影面积，xf、yf、zf分别为水线面漂心的纵向、横向及垂向坐标，i
xf
、i
yf
、i
xyf
分别为水线面对过漂心平行于x、y轴的惯性矩和惯性积，令
[0059]
用长边为船舶总长度l、短边为船舶总宽度b的长方形代替水线面的投影，根据长方形的对称性，则有如下公式：
[0060][0061][0062][0063]
雅克比系数矩阵则简化为：
[0064][0065]
式中，长方形面积s＝lb，
[0066]
则最终简化后雅克比系数矩阵为：
[0067][0068]
d、判断是否满足平衡条件
[0069]
给定船舶初始船中吃水t0、横倾角θ0及纵倾角判断平衡方程组是否满足平衡条件，即判断总载重量和排水量是否相等、重心和浮心是否在同一条垂线上，具体公式如下：
[0070][0071]
如果不满足，则转步骤e；如果满足，则采用纵向梯形积分法计算当前浮态下的船舶排水量v0、浮心纵向坐标x
b0
、浮心横向坐标y
b0
、浮心垂向坐标z
b0
，计算公式为：
[0072][0073]
式中，a0为船长方向x处横剖面的面积，x
f0
、y
f0
、z
f0
分别为船长方向x处横剖面质心纵向坐标、横向坐标及垂向坐标。
[0074]
并按如下公式计算gz0，输出结果，结束。
[0075]
gz0＝(y
b0-yg)cosθ+(z
b0-zg)sinθ
[0076]
e、采用牛顿迭代计算求解平衡方程组
[0077]
e1、进行第一次迭代，计算简化后雅克比系数矩阵；
[0078][0079]
e2、全选主元高斯消去法求解二元一次方程组，得到船中吃水变化量δt0及纵倾角变化量平衡方程组为：
[0080][0081]
即：
[0082]
[0083]
e3、计算当前船舶船中吃水及纵倾角
[0084][0085]
式中，t1为第1次迭代后的船中吃水，为第1次迭代后的纵倾角。
[0086]
e4、采用梯形积分法计算当前浮态下的船舶排水量v1及浮心纵向坐标x
b1
、浮心横向坐标y
b1
及浮心垂向坐标z
b1
。
[0087][0088]
式中，a1为船舶在当前浮态下沿船长方向x处横剖面的面积，x
f1
、y
f1
、z
f1
分别为船长方向x处横剖面质心纵向坐标、横向坐标及垂向坐标。
[0089]
e5、判断是否满足平衡方程组平衡条件精度要求：
[0090][0091]
式中ε1为排水量迭代精度、ε2为纵倾角迭代精度。
[0092]
若满足，按下式计算gz：
[0093]
gz＝(y
b1-yg)cosθ+(z
b1-zg)sinθ
[0094]
输出结果，结束；
[0095]
若不满足，令i＝2，转步骤f继续迭代计算；
[0096]
f、进行第i次迭代求解平衡方程组；
[0097]
f1、采用梯形积分法计算第i次的船舶排水量vi、浮心纵向坐标x
bi
、浮心横向坐标y
bi
、浮心垂向坐标z
bi
；
[0098][0099]
式中，ai为船舶在当前浮态下沿船长方向x处横剖面的面积，x
fi
、y
fi
、z
fi
分别为船
长方向x处横剖面质心纵向坐标、横向坐标及垂向坐标。
[0100]
f2、进行第i次迭代计算雅克比系数矩阵；
[0101][0102]
f3、全选主元高斯消去法求解二元一次方程组，得到第i次的船中吃水变化量δti及纵倾角变化量
[0103][0104]
f4、计算第i次迭代后船中吃水及纵倾角：
[0105][0106]
式中，ti及为第i次船中吃水及船舶纵倾角正切值，t
i-1
及为第i-1次船中吃水及船舶纵倾角正切值。
[0107]
f5、判断是否满足平衡方程组平衡条件精度要求：
[0108][0109]
若满足，按下式计算gz；
[0110]
gz＝(y
bi-yg)cosθ+(z
bi-zg)sinθ
[0111]
并输出当结果，结束；
[0112]
若不满足，令i＝i+1，转步骤f继续迭代计算。
[0113]
与现有技术相比，本发明具有以下有益效果：
[0114]
1、相比于“离线插值法”依赖于船舶装载手册，局限性较大，只能计算船舶纵倾角较小的时候，当吃水差超过装载手册给定范围的时候，计算误差会增大。本发明以船舶3d数据为基础直接计算，不依赖于装载手册，适用于船舶任意浮态，普适性较好。
[0115]
2、本发明相比于“实时计算法”，用长方形代替船舶水线面，减少雅克比系数矩阵的计算量，只需计算出船舶任意浮态下的排水体积及浮心，程序健壮稳定，编程实现简单。
[0116]
3、相比于“优化法”，本发明收敛速度更快，工程应用价值明显。
附图说明
[0117]
图1为gz曲线示意图。
[0118]
图2为本发明的流程图。
[0119]
图3为船体坐标系横向示意图。
[0120]
图4为船体坐标系纵向示意图。
[0121]
图5为水线面投影长方形包围盒。
[0122]
图6为液货船“scaquarius”船舶外壳3d模型示意图。
[0123]
图7为sd18s破损工况示意图。
[0124]
图8为图7的俯视图。
[0125]
图9为破损工况sd18s gz计算结果。
[0126]
图10为sd 55p破损工况示意图。
[0127]
图11为图10的俯视图。
[0128]
图12为破损工况sd 55p gz曲线计算结果。
具体实施方式
[0129]
下面结合附图对本发明做进一步详细地描述，图1所示为gz曲线计算原理示意图、图2为计算流程图、图3和图4为船体坐标系示意图、图5为水线面投影长方形包围盒示意图。本发明以38000dwt液货船“scaquarius”为例进行了实例计算，图6为船舶外壳的3d模型。选取2个典型的破舱工况进行计算，如表1、图7、8、10、11所示。破损工况sd18s的破损舱室包括no.1b.w.t(s)、no.2b.w.t(s)、no.2cargo tank(s)、no.3cargo tank(s)、no.4cargo tank(s)，破损工况sd55p的破损舱室包括rinse fresh water tank(p)、no.13cargo tank(p)、no.2hfo tk(p)、no.1l.s mgo stor tk、pump room、engine room、steer gear room、no.12cargo tank(p)、no.6b.w.t(p)。计算结果见表2、表3、图9、图12，从计算结果可以得到如下结果：本发明所提出的计算方法在计算液货船破损gz曲线时计算精度较高，破损工况sd 55p的各个横倾角度最大计算误差0.03m，平均计算误差0.023m，破损工况sd 18s的各个横倾角度最大计算误差0.02m，平均计算误差0.012m，满足船级社0.5m的精度要求。
[0130]
表1破损组合
[0131]
[0132][0133]
表2破损工况sd18s gz计算误差
[0134]
[0135][0136]
表3破损工况sd 55p gz计算误差
[0137]
横倾角装载手册本发明误差(m)0-0.63-0.6210.0091-0.47-0.4670.00350.150.1550.005100.920.9360.016202.52.5140.014303.823.8370.017404.34.320.02504.084.0950.015603.463.4720.012752.142.1520.012
[0138]
本发明不局限于本实施例，任何在本发明披露的技术范围内的等同构思或者改变，均列为本发明的保护范围。