CN114735152B - 一种基于面元法的全附体舰船不沉性计算方法 - Google Patents

Abstract

本发明涉及一种基于面元法的全附体舰船不沉性计算方法，首先，根据舰船和破损舱室的几何模型，建立相应的三角形面元模型；接着，设置舰船面元模型的法线方向朝外，破损舱室面元模型的法线方向朝内，并以船体坐标系原点为基准将两模型的自由度耦合；然后，输入舰船破损前的重量重心以及初始浮态，基于面元积分计算舰船的排水量、水线面面积等要素；最后，基于牛顿迭代法计算舰船在破损进水后的浮态，提取舰船横倾角、初稳性高、最小稳性力臂等值；在每个迭代步采用面元切分的方法匹配舰船的湿表面面元。本发明能够针对任意船型的全附体舰船开展不沉性计算校核，解决了传统不沉性计算方法在附体、船型方面适用性差的问题，具有广泛的适用性。

Description

一种基于面元法的全附体舰船不沉性计算方法

技术领域

本发明属于舰船安全性技术领域，具体涉及一种基于面元法的全附体舰船不沉性计算方法。

背景技术

舰船在海上执行任务时，在复杂海况和敌方武器的攻击下，易发生破舱进水严重威胁舰船的生命力。因此，舰船需要具备在一定数目水密舱段破损后，仍能保持不沉并保留部分稳性的能力。在舰船设计阶段，需要开展舰船不沉性计算分析，使其满足规范要求。

舰船不沉性计算的核心是求解破舱条件下的舰船浮态，破损舱室通常采用浮力损失法进行处理，即认为破舱不提供浮力；浮态计算通常首先采用二维切线法求解包括排水量、静水面面积等排水要素，在此基础上基于牛顿迭代法进行浮态计算。

传统计算方法对于舰船不沉性问题具有较好的适用性，但同时也存在一定的不足：传统的舰船不沉性计算方法基于二维切线法进行船体建模，由于二维切线对复杂几何体(例如轴支架、螺旋桨等)的拟合程度有限，因此难以精确表达形状复杂的舰船附体；且二维切线法通常采用沿船长方向的线性积分，该方法针对表面光顺的主船体计算精度较高，但是对于存在形状突变(例如舷台结构)的复杂船型，传统计算方法的排水量求解精度具有一定的局限性，因此适用性较差。

发明内容

本发明要解决的技术问题在于针对上述传统不沉性计算方法在附体、船型方面适用性差的问题，提供一种基于面元法的全附体舰船不沉性计算方法，该方法能够准确计算全附体舰船在破舱进水后的浮态，评估舰船的不沉性，且针对任何船型均具有良好的适用性。

本发明为解决上述提出的技术问题所采用的技术方案为：

一种基于面元法的全附体舰船不沉性计算方法，包括以下步骤：

S1、建立面元模型：基于舰船三维模型，采用成熟的面元划分商业软件，建立舰船主船体包括附体的面元模型，面元采用三角形面元，法线方向指向舰船外部；同时建立破损舱室的面元模型，面元采用三角形面元，法线方向指向舰船内部；

S2、指定舰船破损前的重量M，以及其重心在船体坐标系下的坐标值(x_g’,y_g’,z_g’)；所述船体坐标系为坐标原点与船体相连的坐标系，用O’-X’Y’Z’表示，原点O’为舰船基平面、中横剖面与中纵剖面的交点，X’轴沿船长方向，Y’轴沿船宽方向，Z’轴沿吃水方向；

S3、指定舰船初始浮态，通过船舯吃水T₀、横倾角φ₀和纵倾角θ₀来表示；

S4、主船体、破损舱室湿表面面元匹配：

S4.1、将船体坐标系下的面元节点坐标转换到固定坐标系中；所述固定坐标系为原点固定在静水面的坐标系，用O-XYZ表示，静水面为Z＝0的平面；当舰船船舯吃水T为0，横倾角φ为0，纵倾角θ为0时，固定坐标系与船体坐标系重合；

S4.2、遍历主船体、破损舱室的所有面元，提取每个面元内的3个节点的Z坐标值，根据Z坐标值将面元集合分为3类：干表面面元——3个节点的Z坐标值均大于等于0；湿表面面元——3个节点的Z坐标值均小于等于0；跨水线面元——3个节点的Z坐标值部分大于0、部分小于0；

S4.3、以静水面为基平面对跨水线面元进行分割，得到新的湿表面面元，新的湿表面面元与步骤S4.2中的湿表面面元合并构成主船体、破损舱室湿表面的面元集合；

S5、基于面元法计算排水体积要素：

将舰船面元模型和破损舱室面元模型合并为一个计算域从而对破损后的舰船浮态进行求解，破损舱室处的浮力相互抵消从而产生浮力损失的效果；所述排水体积要素包括：排水体积V、排水体积沿X轴的一次积分L_x、排水体积沿Y轴的一次积分L_y、排水体积沿Z轴的一次积分L_z；根据体积分的高斯公式，将排水体积分转化为湿表面的面积分进行求解；

S6、基于面元法计算水线面面积要素：

所述水线面面积要素包括：水线面面积A、水线面面积沿X轴的一次积分S_x、水线面面积沿Y轴的一次积分S_y、水线面面积沿X轴的二次积分I_xx、水线面面积沿Y轴的二次积分I_yy、水线面面积沿X、Y轴的二次积分I_xy；

S7、采用牛顿迭代法计算新的浮态；

S8、判断是否达到平衡状态，若新浮态下舰船满足平衡条件，则进入步骤9；若新浮态下舰船不满足平衡条件，则将计算得到的新浮态作为初始浮态，重新回到步骤S3进行迭代计算，直到满足平衡条件；

S9、根据舰船在破损状态下的平衡浮态，采用舰船静水力稳性计算方法，计算舰船的初稳性高、最小稳性力臂，将计算结果按舰船不沉性规范要求进行校核。

上述方案中，步骤S4.1中，节点坐标转换公式为：

式中，x，y，z表示节点在固定坐标系下的坐标；x’，y’，z’表示节点在船体坐标系下的坐标；T、φ、θ表示舰船的船舯吃水、横倾角和纵倾角。

上述方案中，步骤S4.3中，对跨水线面元进行分割的方法如下：

(1)提取面元的三个节点N1、N2、N3，其Z坐标的乘积为m；

(2)若m＞0，则：以Z＞0的节点为N1，不改变法线方向，调整节点顺序；以N1、N2做直线，与静水面求交点生成N4；以N1、N3做直线，与静水面求交点生成N5；生成新湿表面面元[N2,N5,N4]和[N2,N3,N5]；

若m＜0，则：以Z＜0的节点为N1，不改变法线方向，调整节点顺序；以N1、N2做直线，与静水面求交点生成N4；以N1、N3做直线，与静水面求交点生成N5；生成新湿表面面元[N1,N4,N5]；

若m＝0，则：以Z＝0的节点为N1，不改变法线方向，调整节点顺序；以N2、N3做直线，与静水面求交点生成N4；若N2的Z坐标大于0，则生成新湿表面面元[N1,N4,N3]，若N2的Z坐标小于0，则生成新湿表面面元[N1,N2,N4]。

上述方案中，步骤S5中，排水体积要素的计算式为：

式中，N为主船体和破损舱室湿表面面元数量；S_i第i个面元；x，y，z表示节点在固定坐标系下的坐标。

上述方案中，步骤S6中，静水面面积要素的计算式为：

式中，N为主船体和破损舱室湿表面面元数量；S_i’是第i个面元在静水面的投影；x，y表示节点在固定坐标系下的纵向和横向坐标。

上述方案中，步骤S7中，舰船破损进水后，在固定坐标系内，舰船浮态方程组为：

M-ρV＝0

My_g-ρL_y＝0

Mx_g-ρL_x＝0

式中：M为舰船破损前的重量；ρ为水的密度；V为排水体积；x_g、y_g为舰船重心在固定坐标系下的纵向和横向坐标；L_x、L_y为排水体积沿X轴、Y轴的一次积分；

已知初始船舯吃水T₀、横倾角φ₀和纵倾角θ₀，采用牛顿迭代法计算新的浮态，递推公式为：

式中，T₁、φ₁、θ₁为舰船在新浮态下的船舯吃水、横倾角、纵倾角；J为浮态方程组的雅克比矩阵；V₀为舰船在初始浮态下对应的排水体积；x_g0、y_g0为初始浮态下，舰船重心在固定坐标系下的纵向和横向坐标；L_x0、L_y0为初始浮态下，舰船排水体积沿X轴、Y轴的一次积分；

在初始浮态下，浮态方程组的雅克比矩阵的求解式为：

式中，A₀为舰船在初始浮态下对应的水线面面积；L_z0为初始浮态下，舰船排水体积沿Z轴的一次积分；z_g0为初始浮态下，舰船重心在固定坐标系下的竖向坐标；S_x0、S_y0为初始浮态下，舰船水线面面积沿X轴、Y轴的一次积分；I_xx0、I_yy0为初始浮态下，舰船水线面面积沿X轴、Y轴的二次积分；I_xy0为初始浮态下，舰船水线面面积沿X轴、Y轴的二次积分。

上述方案中，步骤S8中，判断舰船是否达到平衡状态的方法为：计算新浮态下的排水体积V₁、排水体积沿X轴、Y轴的一次积分L_x1、L_y1，设置浮态平衡阈值ε₁、ε₂、ε₃，当达到平衡浮态时，新浮态下舰船重量和排水量的差值的绝对值小于ε₁，舰船重心X轴坐标x_g1和浮心X轴坐标的差值的绝对值小于ε₂，舰船重心Y轴坐标y_g1和浮心Y轴坐标的差值的绝对值小于ε₃；如下式所示：

|M-ρV₁|＜ε₁

式中：M为舰船破损前的重量；ρ为水的密度。

本发明的有益效果在于：

1、本发明建立了一种基于面元法的舰船不沉性计算方法，该方法首先根据舰船和破损舱室的几何模型，建立相应的三角形面元模型；接着设置舰船面元模型的法线方向朝外，破损舱室面元模型的法线方向朝内，并以船体坐标系原点为基准将两模型的自由度耦合；然后输入舰船破损前的重量重心以及初始浮态，基于面元积分计算舰船的排水量、水线面面积等要素；最后基于牛顿迭代法计算舰船在破损进水后的浮态，提取舰船横倾角、初稳性高、最小稳性力臂等值。本发明方法在每个迭代步采用面元切分的方法匹配舰船的湿表面面元，能够针对任意船型的全附体舰船开展不沉性计算校核。

2、本方法解决了传统不沉性计算方法在附体、船型方面适用性差的问题，具有广泛的适用性。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明方法中船体坐标系与固定坐标系的示意图；

图2是本发明方法中三角形面元的静水面分割情况；

图3是本发明方法中跨水线面元分割算法的流程图；

图4是本发明实施例中ONRT标准内倾船三维模型；

图5是本发明实施例中ONRT标准内倾船中纵剖面及液体舱布置图；

图6是本发明实施例中4～6号水密舱段右舷不对称进水的三维破舱模型；

图7是本发明实施例中主船体面元模型；

图8是图7所示主船体面元模型的尾部结构局部放大图；

图9是图8所示尾部结构中螺旋桨和轴支架的局部放大图；

图10是本发明实施例中4～6号水密舱段右舷不对称进水的破舱面元模型；

图11是本发明实施例中ONRT标准内倾船三维模型(不考虑附体)；

图12是本发明实施例中1～4号水密舱段右舷不对称进水的三维破舱模型；

图13是本发明实施例中5～8号水密舱段右舷不对称进水的三维破舱模型；

图14是本发明实施例中9～12号水密舱段右舷不对称进水的三维破舱模型；

图15是本发明实施例中13～14号水密舱段右舷不对称进水的三维破舱模型。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

本发明提供了一种基于面元法的全附体舰船不沉性计算方法，包括以下步骤：

S1、建立面元模型。基于舰船三维模型，采用成熟的面元划分商业软件，建立舰船主船体包括附体的面元模型，面元采用三角形面元，法线方向指向舰船外部；同时建立破损舱室的面元模型，面元采用三角形面元，法线方向指向舰船内部。

S2、指定舰船破损前的重量M，以及其重心在船体坐标系下的坐标值(x_g’,y_g’,z_g’)。所述船体坐标系为坐标原点与船体相连的坐标系，用O’-X’Y’Z’表示，原点O’为舰船基平面、中横剖面与中纵剖面的交点，X’轴沿船长方向，Y’轴沿船宽方向，Z’轴沿吃水方向，参见图1。

S3、指定舰船初始浮态。可通过船舯吃水T₀、横倾角φ₀和纵倾角θ₀来表示。

S4、主船体、破损舱室湿表面面元匹配。舰船浮力和破损舱室损失的浮力仅与其湿表面有关，随着舰船浮态的变化，其湿表面也会产生相应的变化，即湿表面面元为求解域。因此需要开展主船体、破损舱室湿表面面元匹配。任意三角形面元与静水面有且仅有三种位置关系：(1)面元在静水面以上；(2)面元在静水面以下；(3)面元跨静水面，一部分在静水面以上，一部分在静水面以下。因此，湿表面面元匹配的流程为：

S4.1、将船体坐标系下的面元节点坐标转换到固定坐标系中。所述固定坐标系为原点固定在静水面的坐标系，用O-XYZ表示，继续参见图1，在固定坐标系，静水面为Z＝0的平面，有利于开展湿表面面元匹配的操作，当舰船船舯吃水T为0，横倾角φ为0，纵倾角θ为0时，固定坐标系与船体坐标系重合。节点坐标转换公式为：

式中，x，y，z表示节点在固定坐标系下的坐标；x’，y’，z’表示节点在船体坐标系下的坐标；T、φ、θ表示舰船的船舯吃水、横倾角和纵倾角。

S4.2、遍历主船体、破损舱室的所有面元，提取每个面元内的3个节点的Z坐标值，根据Z坐标值将面元集合分为3类：干表面面元——3个节点的Z坐标值均大于等于0；湿表面面元——3个节点的Z坐标值均小于等于0；跨水线面元——3个节点的Z坐标值部分大于0、部分小于0。

S4.3、以静水面(Z＝0)为基平面对跨水线面元进行分割，考虑三角形面元与静水面的位置关系，存在3种分割情况，如图2所示：

(1)1个节点在静水面上方，2个节点在静水面下方；

(2)2个节点在静水面上方，1个节点在静水面下方；

(3)1个节点在静水面上方，1个节点在静水面下方，1个节点在静水面内。

针对不同的面元分割类型，分别求解面元边线与静水面的交点生成1～2个新增节点，在此基础上，按照原始面元节点的连接方向将其分割为2～3个子面元，确保子面元的法线方向不变，从而分割提取出新的湿表面面元。其计算的具体步骤如图3所示，对跨水线面元进行分割的方法如下：

(1)提取面元的三个节点N1、N2、N3，其Z坐标的乘积为m；

(2)若m＞0，则：以Z＞0的节点为N1，不改变法线方向，调整节点顺序；以N1、N2做直线，与静水面求交点生成N4；以N1、N3做直线，与静水面求交点生成N5；生成新湿表面面元[N2,N5,N4]和[N2,N3,N5]；

若m＜0，则：以Z＜0的节点为N1，不改变法线方向，调整节点顺序；以N1、N2做直线，与静水面求交点生成N4；以N1、N3做直线，与静水面求交点生成N5；生成新湿表面面元[N1,N4,N5]；

若m＝0，则：以Z＝0的节点为N1，不改变法线方向，调整节点顺序；以N2、N3做直线，与静水面求交点生成N4；若N2的Z坐标大于0，则生成新湿表面面元[N1,N4,N3]，若N2的Z坐标小于0，则生成新湿表面面元[N1,N2,N4]。

得到的新的湿表面面元与步骤S4.2中的湿表面面元合并构成主船体、破损舱室湿表面的面元集合。

S5、基于面元法计算排水体积要素。在固定坐标系下，因为破损舱室的面元模型法线方向指向船内，其计算得到的“浮力”方向向下，因此，可以将舰船面元模型和破损舱室面元模型合并为一个计算域从而对破损后的舰船浮态进行求解，破损舱室处的浮力相互抵消从而产生浮力损失的效果。采用牛顿迭代法求解舰船浮态需要预先计算4个排水体积要素，包括：排水体积V、排水体积沿X轴的一次积分L_x、排水体积沿Y轴的一次积分L_y、排水体积沿Z轴的一次积分L_z。根据体积分的高斯公式，可以将排水体积分转化为湿表面的面积分进行求解，推导排水体积要素的计算式为：

式中，N为主船体和破损舱室湿表面面元数量；S_i第i个面元；x，y，z表示节点在固定坐标系下的坐标。

S6、基于面元法计算水线面面积要素。采用牛顿迭代法求解舰船浮态需要预先计算6个水线面面积要素，包括：水线面面积A、水线面面积沿X轴的一次积分S_x、水线面面积沿Y轴的一次积分S_y、水线面面积沿X轴的二次积分I_xx、水线面面积沿Y轴的二次积分I_yy、水线面面积沿X、Y轴的二次积分I_xy。静水面面积要素的计算式为：

式中，N为主船体和破损舱室湿表面面元数量；S_i’是第i个面元在静水面的投影；x，y表示节点在固定坐标系下的纵向和横向坐标。

S7、采用牛顿迭代法计算新的浮态。舰船破舱进水后，当其在固定坐标系下保持平衡时，舰船重力和浮力大小相等方向相反，且作用在同一条直线上。在固定坐标系内，舰船浮态方程组为：

M-ρV＝0

My_g-ρL_y＝0

Mx_g-ρL_x＝0

式中：M为舰船破损前的重量；ρ为水的密度；V为排水体积；x_g、y_g为舰船重心在固定坐标系下的纵向和横向坐标；L_x、L_y为排水体积沿X轴、Y轴的一次积分。

已知初始船舯吃水T₀、横倾角φ₀和纵倾角θ₀，采用牛顿迭代法计算新的浮态，递推公式为：

式中，T₁、φ₁、θ₁为舰船在新浮态下的船舯吃水、横倾角、纵倾角；J为浮态方程组的雅克比矩阵；V₀为舰船在初始浮态下对应的排水体积；x_g0、y_g0为初始浮态下，舰船重心在固定坐标系下的纵向和横向坐标；L_x0、L_y0为初始浮态下，舰船排水体积沿X轴、Y轴的一次积分。

在初始浮态下，浮态方程组的雅克比矩阵的求解式为：

式中，A₀为舰船在初始浮态下对应的水线面面积；L_z0为初始浮态下，舰船排水体积沿Z轴的一次积分；z_g0为初始浮态下，舰船重心在固定坐标系下的竖向坐标；S_x0、S_y0为初始浮态下，舰船水线面面积沿X轴、Y轴的一次积分；I_xx0、I_yy0为初始浮态下，舰船水线面面积沿X轴、Y轴的二次积分；I_xy0为初始浮态下，舰船水线面面积沿X轴、Y轴的二次积分。

S8、判断是否达到平衡状态。计算新浮态下的排水体积V₁、排水体积沿X轴、Y轴的一次积分L_x1、L_y1，设置浮态平衡阈值ε₁、ε₂、ε₃，当达到平衡浮态时，新浮态下舰船重量和排水量的差值的绝对值小于ε₁，舰船重心X轴坐标x_g1和浮心X轴坐标的差值的绝对值小于ε₂，舰船重心Y轴坐标y_g1和浮心Y轴坐标的差值的绝对值小于ε₃。如下式所示：

|M-ρV₁|＜ε₁

若新浮态下舰船满足平衡条件，则进入步骤9；若新浮态下舰船不满足平衡条件，则将计算得到的新浮态作为初始浮态，重新回到步骤S3进行迭代计算，直到满足平衡条件。

S9、根据舰船在破损状态下的平衡浮态，采用舰船静水力稳性计算方法，计算舰船的初稳性高、最小稳性力臂，将计算结果按舰船不沉性规范要求进行校核。

以下以ONRT标准内倾船型为例，对本发明方法进行说明，该实施例不构成对本发明的限制。内倾船的三维几何模型如图4所示，由主船体和附体组成。附体包括舭龙骨、呆木、轴包套、尾轴、轴支架、螺旋桨、舵，其相应位置已在图4中标出。该船的设计吃水5.55m，设计水线长154m，型深14.5m。

沿船长方向均匀设置13道水密舱壁，将该舰船平均分隔为14个水密舱段，沿船首至船尾方向依次编号为1～14号水密舱段，2～13号水密舱段底部设置有液体舱，考虑到舰船艏部狭窄，在艏部2～3号水密舱段底部各布置1个液体舱，4～13号水密舱段底部各对称布置2个液体舱。距基线4.5m高度做剖面1以方便显示液体舱布置。中纵剖面及液体舱布置如图5所示。

基于本发明提出的基于面元法的全附体舰船不沉性计算方法，开展舰船4～6号水密舱段同时发生右舷不对称进水情况下的不沉性计算。4～6号水密舱段右舷不对称进水时的破损舱的三维模型如图6所示。

S1.建立面元模型。采用成熟的商业软件，分别建立舰船主船体的面元模型和破损水密舱段的面元模型，采用三角形面元，主船体的面元法线方向指向舰船外部，破损水密舱段的面元法线方向指向舰船内部。针对螺旋桨、轴支架等复杂曲面，加密面元数量使其良好贴合曲面，提高计算精度。主船体的面元数量为103568，4号破损水密舱段的面元数量为848，5号破损水密舱段的面元数量为1112，6号破损水密舱段的面元数量为1094。主船体的面元模型如图7-9所示，破损水密舱段的面元模型如图10所示。

S2.设定本舰在未破损状态下的重量为7500t，在船体坐标系下的重心坐标为(-2.40,0.00,7.55)，海水密度设置为1.025t/m³。

S3.指定初始浮态。设定本舰的初始浮态为：船舯吃水T₀＝5.00m，横倾角φ₀＝0°，纵倾角θ₀＝0°。

S4.主船体、破损舱室湿表面面元匹配。主船体和破损水密舱段总的面元数量为106622，在初始浮态下，通过湿表面面元匹配算法计算得到主船体、破损水面舱段的湿表面面元集合，面元数量为94465。

S5.基于面元法计算排水体积要素。采用面元法计算得到的4个排水体积要素为：排水体积V＝5340.58m³、排水体积沿X轴的一次积分L_x＝-59711.74m⁴、排水体积沿Y轴的一次积分L_y＝0.00m⁴、排水体积沿Z轴的一次积分L_z＝-10870.48m⁴。

S6.基于面元法计算水线面面积要素。采用面元法计算得到的6个水线面面积要素为：水线面面积A＝1622.70m²、水线面面积沿X轴的一次积分S_x＝-33892.74m³、水线面面积沿Y轴的一次积分S_y＝0.00m³、水线面面积沿X轴的二次积分I_xx＝2561304.38m⁴、水线面面积沿Y轴的二次积分I_yy＝38079.21m⁴、水线面面积沿X、Y轴的二次积分I_xy＝1.47m⁴。

S7.采用牛顿迭代法计算新的浮态。采用牛顿迭代法计算新的浮态为：船舯吃水T₁＝7.18m，横倾角φ₁＝-14.57°，纵倾角θ₀＝2.64°。

S8.判断是否达到平衡状态。设置浮态平衡阈值ε₁＝1t、ε₂＝0.01m、ε₃＝0.01m，新浮态对应的排水量为5474t，不满足平衡条件，将计算得到的浮态设置为初始浮态，重新回到步骤S3进行迭代计算。迭代计算6次后，达到平衡状态，迭代计算过程如下表1所示。

表1基于牛顿迭代法的破损舰船浮态过程

S9.基于牛顿迭代法计算得到的舰船破损状态下的浮态为：船舯吃水T₁＝7.29m，横倾角φ₁＝2.46°，纵倾角θ₀＝3.03°。根据舰船在破损状态下的平衡浮态，采用舰船静水力稳性计算方法，计算得到舰船的初稳性高为1.11m，最小初稳性臂为0.64m，将计算结果按舰船不沉性规范要求进行校核，满足规范要求。

为了验证本发明方法的准确性，下面将本计算方法与成熟商业软件MAXSURF进行对比计算验证。分别计算ONRT标准内倾船的1～14号水密舱段发生右舷不对称进水后的浮态，设定本舰在未破损状态下的重量为7500t，在船体坐标系下的重心坐标为(-2.00,0.00,7.55)，海水密度设置为1.025t/m³。考虑到MAXSURF的浮态计算通常仅适用于裸船体，因此，准确性验证计算未将附体纳入计算域。主船体的三维模型如图11所示，14个破损舱段的三维模型如图12～图15所示。

舰船不同水密舱段破损进水后的浮态计算对比如表2所示，由表2可知，该方法的结算结果与MAXSURF基本一致，船舯吃水最大误差0.01m，横倾角的最大误差为0.17°，纵倾角的最大误差为0.01°。表明本发明计算方法具有较好的准确性。

表2本方法与MAXSURF计算结果对比

上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (6)

1.一种基于面元法的全附体舰船不沉性计算方法，其特征在于，包括以下步骤：

S1、建立面元模型：基于舰船三维模型，采用成熟的面元划分商业软件，建立舰船主船体包括附体的面元模型，面元采用三角形面元，法线方向指向舰船外部；同时建立破损舱室的面元模型，面元采用三角形面元，法线方向指向舰船内部；

S2、指定舰船破损前的重量M，以及其重心在船体坐标系下的坐标值（x_g’,y_g’,z_g’）；所述船体坐标系为坐标原点与船体相连的坐标系，用O’-X’Y’Z’表示，原点O’为舰船基平面、中横剖面与中纵剖面的交点，X’轴沿船长方向，Y’轴沿船宽方向，Z’轴沿吃水方向；

S3、指定舰船初始浮态，通过船舯吃水T₀、横倾角φ₀和纵倾角θ₀来表示；

S4、主船体、破损舱室湿表面面元匹配：

S4.1、将船体坐标系下的面元节点坐标转换到固定坐标系中；所述固定坐标系为原点固定在静水面的坐标系，用O-XYZ表示，静水面为Z=0的平面；当舰船船舯吃水T为0，横倾角φ为0，纵倾角θ为0时，固定坐标系与船体坐标系重合；

S4.2、遍历主船体、破损舱室的所有面元，提取每个面元内的3个节点的Z坐标值，根据Z坐标值将面元集合分为3类：干表面面元——3个节点的Z坐标值均大于等于0；湿表面面元——3个节点的Z坐标值均小于等于0；跨水线面元——3个节点的Z坐标值部分大于0、部分小于0；

S4.3、以静水面为基平面对跨水线面元进行分割，得到新的湿表面面元，新的湿表面面元与步骤S4.2中的湿表面面元合并构成主船体、破损舱室湿表面的面元集合；对跨水线面元进行分割的方法如下：

（1）提取面元的三个节点N1、N2、N3，其Z坐标的乘积为m；

（2）若m＞0，则：以Z＞0的节点为N1，不改变法线方向，调整节点顺序；以N1、N2做直线，与静水面求交点生成N4；以N1、N3做直线，与静水面求交点生成N5；生成新湿表面面元[N2,N5,N4]和[N2,N3,N5]；

若m＜0，则：以Z＜0的节点为N1，不改变法线方向，调整节点顺序；以N1、N2做直线，与静水面求交点生成N4；以N1、N3做直线，与静水面求交点生成N5；生成新湿表面面元[N1,N4,N5]；

若m=0，则：以Z=0的节点为N1，不改变法线方向，调整节点顺序；以N2、N3做直线，与静水面求交点生成N4；若N2的Z坐标大于0，则生成新湿表面面元[N1,N4,N3]，若N2的Z坐标小于0，则生成新湿表面面元[N1,N2,N4]；

S5、基于面元法计算排水体积要素：

将舰船面元模型和破损舱室面元模型合并为一个计算域从而对破损后的舰船浮态进行求解，破损舱室处的浮力相互抵消从而产生浮力损失的效果；所述排水体积要素包括：排水体积V、排水体积沿X轴的一次积分L_x、排水体积沿Y轴的一次积分L_y、排水体积沿Z轴的一次积分L_z；根据体积分的高斯公式，将排水体积分转化为湿表面的面积分进行求解；

S6、基于面元法计算水线面面积要素：

所述水线面面积要素包括：水线面面积A、水线面面积沿X轴的一次积分S_x、水线面面积沿Y轴的一次积分S_y、水线面面积沿X轴的二次积分I_xx、水线面面积沿Y轴的二次积分I_yy、水线面面积沿X、Y轴的二次积分I_xy；

S7、采用牛顿迭代法计算新的浮态；

S8、判断是否达到平衡状态，若新浮态下舰船满足平衡条件，则进入步骤9；若新浮态下舰船不满足平衡条件，则将计算得到的新浮态作为初始浮态，重新回到步骤S3进行迭代计算，直到满足平衡条件；

S9、根据舰船在破损状态下的平衡浮态，采用舰船静水力稳性计算方法，计算舰船的初稳性高、最小稳性力臂，将计算结果按舰船不沉性规范要求进行校核。

2.根据权利要求1所述的基于面元法的全附体舰船不沉性计算方法，其特征在于，步骤S4.1中，节点坐标转换公式为：

式中，x，y，z表示节点在固定坐标系下的坐标；x’，y’，z’表示节点在船体坐标系下的坐标；T、φ、θ表示舰船的船舯吃水、横倾角和纵倾角。

3.根据权利要求1所述的基于面元法的全附体舰船不沉性计算方法，其特征在于，步骤S5中，排水体积要素的计算式为：

式中，N为主船体和破损舱室湿表面面元数量；S_i第i个面元；x，y，z表示节点在固定坐标系下的坐标。

4.根据权利要求1所述的基于面元法的全附体舰船不沉性计算方法，其特征在于，步骤S6中，静水面面积要素的计算式为：

式中，N为主船体和破损舱室湿表面面元数量；S_i’是第i个面元在静水面的投影；x，y表示节点在固定坐标系下的纵向和横向坐标。

5.根据权利要求1所述的基于面元法的全附体舰船不沉性计算方法，其特征在于，步骤S7中，舰船破损进水后，在固定坐标系内，舰船浮态方程组为：

式中：M为舰船破损前的重量；ρ为水的密度；V为排水体积； x_g、y_g为舰船重心在固定坐标系下的纵向和横向坐标；L_x、L_y为排水体积沿X轴、Y轴的一次积分；

已知初始船舯吃水T₀、横倾角φ₀和纵倾角θ₀，采用牛顿迭代法计算新的浮态，递推公式为：

式中，T₁、φ₁、θ₁为舰船在新浮态下的船舯吃水、横倾角、纵倾角；J为浮态方程组的雅克比矩阵；V₀为舰船在初始浮态下对应的排水体积；x_g0、y_g0为初始浮态下，舰船重心在固定坐标系下的纵向和横向坐标；L_x0、L_y0为初始浮态下，舰船排水体积沿X轴、Y轴的一次积分；

在初始浮态下，浮态方程组的雅克比矩阵的求解式为：

式中，A₀为舰船在初始浮态下对应的水线面面积；L_z0为初始浮态下，舰船排水体积沿Z轴的一次积分；z_g0为初始浮态下，舰船重心在固定坐标系下的竖向坐标；S_x0、S_y0为初始浮态下，舰船水线面面积沿X轴、Y轴的一次积分；I_xx0、I_yy0为初始浮态下，舰船水线面面积沿X轴、Y轴的二次积分；I_xy0为初始浮态下，舰船水线面面积沿X轴、Y轴的二次积分。

6.根据权利要求1所述的基于面元法的全附体舰船不沉性计算方法，其特征在于，步骤S8中，判断舰船是否达到平衡状态的方法为：计算新浮态下的排水体积V₁、排水体积沿X轴、Y轴的一次积分L_x1、L_y1，设置浮态平衡阈值ε₁、ε₂、ε₃，当达到平衡浮态时，新浮态下舰船重量和排水量的差值的绝对值小于ε₁，舰船重心X轴坐标x_g1和浮心X轴坐标的差值的绝对值小于ε₂，舰船重心Y轴坐标y_g1和浮心Y轴坐标的差值的绝对值小于ε₃；如下式所示：

式中：M为舰船破损前的重量；ρ为水的密度。