船舶大幅横摇运动的三维数值模拟方法与流程

本发明涉及的是一种模拟方法，特别涉及一种利用三维泰勒展开边界元方法模拟船舶大幅横摇运动的方法。

背景技术：

首尾变化剧烈的船舶，在特定的海浪工况下容易发生参数横摇现象，即大幅横摇运动。特别是集装箱船极容易发生该现象。船舶一旦发生参数横摇现象则会对货物带来损失，对船舶安全带来影响。因此船舶参数横摇现象发生概率统计是水动力学研究的热点。

船舶运动预报数值模拟中需要考虑定常叠模势对非定常速度势的影响。因此自由面条件实施对数值预报精度也有影响。目前常用的自由面条件包括nk自由面条件和db自由面条件。因为db自由面条件设计叠模势二阶导数的计算，因此要比nk自由面条件要复杂。而尖角边界处速度势二阶导数精确求解是数值模拟的难点。

集装箱船舶尾部均为大外飘船型，利用时域格林函数方法求解存在数值发散问题。因此对于大外飘船型运动时域预报均采用简单格林函数法。对于远方辐射条件实施有阻尼区法和时域格林函数法。相比于阻尼区法，时域格林函数方法大大减少自由面网格数目，能够提高数值模拟的效率。

技术实现要素：

本发明的目的在于提供一种能够为船舶参数横摇现象发生概率统计提供基础的船舶大幅横摇运动的三维数值模拟方法。

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

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

步骤2，计算泰勒展开边界元方法所需的边界积分方程所涉及的影响系数矩阵；

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

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

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

步骤6，不规则波分解，线性叠加成入射波时历；

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

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

步骤9，船舶大幅运动预报方程建模，采用四阶龙格库塔方法步进求解运动方程，进行船舶在顶浪或者斜浪中非线性运动评估；

步骤10，依据船舶大幅运动响应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)将波浪谱进行傅里叶变换，得到对应波浪谱的子波参数；本方法将输入的波浪时历利用ittc双参数谱的波谱公式通过有义波高和周期表示如下：

其中，t1为谱心周期；h1/3为有义波高，ω为圆频率，s(ω)为海浪谱密度。

利用等分能量法，将谱能量分为num份，各组成波的波幅如下：

其中ω1～ω2为截断频率，利用线性叠加原理即可得到入射波浪的时历表达式：

其中εi为各组成波的相位，而且

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

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

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

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

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

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

静水压力为：

ps＝-ρgz

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

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

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

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

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

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

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

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

10)计算船舶大幅运动响应rao，进行不规则波中船舶大幅运动数值模拟和特征统计。

根据船舶运动方程时域模拟，获得船舶运动响应时历，由此利用修正的周期图平均法(也称为welch方法)获得船舶航行中遭遇的运动响应谱si(ω)和海浪谱sζ(ω)

其中l为将时历分的段数，ii(ω)为加了海明窗后每一段的修正周期图。

假设运动响应谱和海浪谱之间满足线性传递函数，按照谱分析基本理论，可获得船舶运动响应rao，如下式所示。

上式中si(ω)代表由某一自由度响应曲线得到的频谱函数，sζ(ω)代表波面的频谱函数。

利用逐波分析理论对模拟的运动响应时历有义值、运动极值和振荡统计进行分析。