预报恶劣海况下低速肥大型船舶波浪增阻的三维数值方法与流程

本发明涉及的是一种数值模拟方法，特别涉及一种利用三维泰勒展开边界元方法模拟低航速肥大型船舶恶劣海况中波浪增阻的数值方法。

背景技术：

船舶波浪增阻值是评估船舶co2排放等级的关键指标，在新船设计开发过程中为满足imo国际海事组织的船舶co2排放指数，以及运力过剩的情况下。船舶设计人员选择降低设计航速，从而降低主机功率。以便满足船舶co2排放指数。对于油船，散货船等肥大型船舶主机储备功率不足，会导致恶劣海况下船舶安全性能降低。因此imo国际海事组织正在讨论制定船舶最小装机功率评估导则。导则核心即如何评估恶劣海况下船舶增阻对航向稳定性的影响。

船舶运动预报数值模拟中需要考虑定常叠模势对非定常速度势的影响。因此自由面条件实施对数值预报精度也有影响。目前常用的自由面条件包括nk自由面条件和db自由面条件。因为db自由面条件设计叠模势二阶导数的计算，因此要比nk自由面条件要复杂。而尖角边界处速度势二阶导数精确求解是数值模拟的难点。imo临时导则规定8级海况。因此需要考虑船舶在恶劣海况中大幅运动对波浪增阻的影响，这是数值预报精度的关键点。

技术实现要素：

本发明的目的在于提供一种能够为船舶最小装机功率评估提供基础的预报恶劣海况下低速肥大型船舶波浪增阻的三维数值方法。

本发明的目的是这样实现的：

步骤1，读取网格文件，利用网格信息进行船舶静水力计算；

步骤2，计算泰勒展开边界元方法所需的边界积分方程所涉及的影响系数矩阵；为了提高数值效率包括关于轴对称面xoz的镜像点源产生的影响系数矩阵。

步骤3，叠模速度势及其空间一阶、二阶导数和mj项求解；

步骤4，时域格林函数及其空间法向导数求解；

步骤5，横摇阻尼系数计算；

步骤6，积分格式自由面条件时域步进自由面离散网格中心点处的速度势；

步骤7，泰勒展开边界元法直接时域扰动波浪力计算；

步骤8，基于弱散射理论入射波浪力(froude-krylov力)与静水回复力计算；

步骤9，船舶大幅运动预报方程建模，采用四阶龙格库塔方法步进求解运动方程；

步骤10，船舶波浪增阻计算；

步骤11，依据船舶运动响应rao和增阻rao，进行不规则波中船舶大幅运动和波浪增阻谱分析计算。

本发明提供了一种预报肥大型船舶(油船、散货船)恶劣海况中波浪增阻数值方法，为船舶最小装机功率评估提供基础。

船舶波浪增阻是二阶波浪载荷，涉及定常和非定常速度势的二阶导数。因此精确求解尖角边界处叠模势，非定常扰动速度势切向诱导速度及高阶导数是有航速定解问题数值模拟的关键。本发明利用泰勒展开边界元方法能够精确求解任意流域边界的二阶导数。弱散射理论可以考虑瞬时波面下入射波浪力和恢复力，考虑船舶大幅运动。

利用本发明提出的方法能够预报集装箱船迎浪规则波中大幅运动rao，波浪增阻rao，以及不规则波中船舶大幅运动和波浪增阻谱分析结果。

附图说明

图1是本发明的流程图。

具体实施方式

下面举例对本发明做更详细的描述。

1)读取网格文件自动提取船舶水线信息(包括首尾驻点空间坐标，水线拟合曲线函数)。结合流场匹配边界智慧参数和水线信息，自动生成符合边界元方法的水面和匹配面离散网格。并基于船体离散网格，计算船舶排水体积，浮心，漂心，惯性矩，湿表面积等静水力参数。检查静水力参数数值结果与物理船舶参数的误差，以此检验船舶网格质量。

2)本发明利用泰勒展开边界元求解各速度势成分及其空间一阶、二阶导数。泰勒展开核心思想是基于格林第三公式形成的边界积分方程进行数值离散求解的方法。对于三维问题，将浮体湿表面离散为若干四边形或三角形单元，在每一单元上，取单元节点坐标均值为中心，在面元中点对偶极强度作泰勒展开并保留一阶导数项，对源强作泰勒展开只保留一阶导数项。并引入场点的切向一阶导数来封闭方程组，从而构成了关于偶极强度、偶强的一阶导数为未知数，源强为已知变量的线性代数方程组。其中偶强切向一阶导数的影响系数包含主值(归一化后为正/负二分之一)。上述操作方法产生的利用边界单元求解边界积分方程数值解的方法称为泰勒展开边界元方法。对于任意面元i可得到如下简化的一阶泰勒展开边界元方法的离散方程组，i＝1,2,…,n，

上式各矩阵中元素表达式：

式中：上角标i和j表示面元编号。以矩阵中某一元素表达式为例做一解释：如：

该方法可同时求解速度势及两个相互正交的切向方向导数。在利用物面法向不可穿透条件，即构成了当地局部坐标系的速度场。可实现速度场在局部坐标下及大地坐标系下的转换。数值结果证明，该方法可明显改善流域边界拐角处的切向诱导速度的计算精度。引入辅助函数再次利用泰勒展开边界元方法求解φ的空间，即速度势空间二阶导数。因此该方法涉及到若干影响系数矩阵计算。

3)总速度势可分解为定常速度势、非定常入射势、辐射速度势和绕射速度势。即：

定常速度势又可分解为来流速度势和叠模速度势，即：

φb＝-ux+φ(3)

非定常势定解问题的物面条件涉及定常势影响，即mj项。对于db假设，其边值定解问题为：

基于db线性假设mj项为：

式中：为物面各点处的位移。为船体平动位移，为船体转动位移，为船体湿表面上各点位置矢径。

4)本发明利用匹配法外传扰动波浪，在匹配面上利用时域格林函数构造匹配条件。时域格林函数为：

5)横摇阻尼计算

利用ikeda方法计算船舶横摇阻尼系数，注意此时得到的总阻尼为无因次等效线性阻尼系数。该方法中，按如下方式无因次化：

其中ρ为水密度，为排水体积，b为型宽，等效线性阻尼系数可表示为：

其中b441为线性阻尼系数，b443为立方阻尼系数，η4为横摇幅值，ωη为横摇固有频率。

为了更准确的求解非线性阻尼系数，本发明分别求出不同横摇幅值下的等效线性阻尼系数，然后根据上式进行最小二乘拟合，即可求得以上非线性阻尼系数。

6)本发明采用时域直接求解非定常扰动速度势，辐射速度势和绕射速度势一起求解。非定常扰动速度势定解问题如式(9)，利用泰勒展开边界元法时域步进求解该初边值问题。

式中自由面条件右端项f的表达式表示如下：

采用积分格式自由面条件步进自由面上各离散单元中心点处的速度势。以任意函数f(t)为例，阐述积分格式自由面条件的核心思想：即对被积函数作时间二次积分。

同理对自由面条件作时间二次积分得：

利用梯形法计算积分格式自由面条件。从而实现自由面条件时间步进。一旦扰动速度势求解完后，利用伯努利方程在平均湿表面上积分，可得到扰动波浪载荷。

7)瞬时波面以下入射波浪力(froude-krylov力)与静水回复力计算；由伯努利方程，可知静水面以下入射波压力为：

静水压力为：

ps＝-ρgz

由于对于入射波力和静水回复力的计算积分至瞬时湿表面，故存在入射波波面在静水面以上的情况，因此本发明采用wheeler法评估静水面以上入射波压力的贡献：

pis＝-ρgz+ρgζ(t)e^k(z-ζ(t))(z≤ζ(t))(16)

此方法保证了波浪压力随深度增加而呈指数衰减。

为了计算瞬时湿表面下静水恢复力和入射波浪力，先将面网格沿船体表面划分至甲板，采用四边形和三角形网格。在运动方程计算得到瞬时船体运动姿态后，通过转换关系将原始正浮状态的网格信息转换为瞬时姿态下的网格信息。在每个时间步判断船体上的湿表面网格在入射波面以下，以上还是与入射波面相截。瞬时网格中，水面以下的部分即为瞬时湿表面。

在获得了瞬时入射波面与船体湿表面交界面以下的瞬时湿表面后，将入射波和静水压力pis在瞬时船体湿表面积分，就获得了船舶受到的非线性入射波力和静水恢复力,表示如下：

在上式中，入射波和静水恢复力矩相对船舶重心的随体坐标系取矩。其原因是后续的船舶运动转动方程是在相对于船舶重心的随体坐标系下建立起来的。

8)根据船舶大幅运动预报方程建模，采用运动方程可进行船舶在顶浪中非线性运动评估。依据牛顿第二定律，船舶垂荡，横摇，纵摇三自由度耦合运动方程为：

利用四阶龙格库塔方法步进求解运动方程。

9)当速度势及其空间一阶、二阶导数计算完成，船舶六自由度运动时历信号计算完成后，利用近场压力积分公式，保留二阶压力载荷项，如下式所示即为船舶波浪增阻时历表达式。对该式取平均值即可得到船舶波浪增阻值。

10)得到船舶六自由度运动rao和波浪增阻rao后，利用谱分析方法能够得到各级海况下船舶运动和波浪增阻的统计值。ittc双参数谱如下式所示：

其中，t1为谱心周期；h1/3为有义波高，ω为圆频率，s(ω)为海浪谱密度。将六自由度运动rao值带入谱分析公式，便得到船舶在实际海况中的运动幅值。