基于泰勒展开边界元法的波浪中船舶回转运动预报方法

1.本发明涉及船舶水动力操纵性领域，特别涉及一种基于泰勒展开边界元法的波浪中船舶回转运动预报方法。

背景技术：

2.近年来，出于对船舶能效指标和航行安全性考虑，船舶耦合水动力性能分析成为研究热点。在船舶水动力学发展初期，受制于理论模型、计算机能力等，不得不进行若干假设、简化，进而开展单一的船舶水动力性能研究。但实际海上航行中，船舶耐波运动和操纵运动是同时发生且相互影响的。ittc委员会专门成立了波浪中船舶操纵委员会，关注船舶操纵与耐波耦合运动分析的前沿进展。
3.目前流行的船舶操纵与耐波运动耦合分析研究方法大致可分为以下四种：船模试验方法；基于粘流cfd的数值模拟方法；基于耐波/操纵统一理论的数值模拟方法；基于耐波/操纵双时间尺度耦合的数值模拟方法。
4.实验方法虽然是最可靠的方法但是消耗的成本较大；利用cfd方法模拟波浪中船舶操纵运动时计算效率较低；而基于统一理论的数值模拟方法虽然引入了二阶波浪力的作用，但相比较双时间尺度方法，二阶波浪力的计算精度较低。
5.双时间尺度模型假设船舶耐波运动是波频运动，可利用三维有航速势流理论方法预报船舶在规则/非规则波中的漂移载荷，而操纵运动是低频运动，可利用mmg方程构建船舶操纵模型从而实现船舶操纵与耐波耦合运动分析模型。在双时间尺度模型中，使用泰勒展开边界元法进行耐波性计算克服了传统常值面元法在求解非光滑边界处切向诱导速度及高阶边界元和其他数值方法在处理拐角处速度势高阶导数精度低的缺点。

技术实现要素：

6.本发明提供了一种基于泰勒展开边界元法的波浪中船舶回转运动预报方法，提高了现有船舶回转运动轨迹精度，实现了船舶长时间数值模拟的稳定性。
7.为实现上述目的，本发明提供的基于泰勒展开边界元法的波浪中船舶回转运动预报方法，包括以下步骤：
8.步骤1:读取网格文件，利用网格信息进行静水力计算；
9.步骤2:通过时域泰勒展开边界元法计算船舶初始状态下的叠模速度势、扰动速度势及其空间一阶、二阶导数；
10.步骤3:直接在船体湿表面压力积分获得一阶力，进而通过船舶6自由度运动方程求得船舶的波频运动，并且通过近场积分方法求得带漂角船舶所受的二阶波浪漂移力和力矩；
11.步骤4：利用步骤3获得的船舶6自由度波频运动更新扰动势边值问题，并且再次利用时域泰勒展开边界元法对扰动势进行求解，获得新时刻的船舶6自由度波频运动，进而求得新时刻带漂角船舶所受的二阶波浪漂移力和力矩。步骤2到步骤4称为船舶耐波性问题，
并且对耐波性问题进行n
t
步求解；
12.步骤5：将获得的二阶波浪漂移力和力矩进行最小二乘拟合，并带入mmg操纵运动模型，计算下一时刻的船舶运动信息；
13.步骤6：根据新的船舶运动信息重新计算叠模速度势和扰动速度势；
14.步骤7:更新船舶的6自由度运动和船舶所受的二阶波浪漂移载荷，从而计算出下一时刻的船舶运动信息。
15.具体而言，所述步骤1具体包括：
16.步骤101:读取网格文件，自动提取船舶水线信息和自动生成符合边界元方法的水面离散网格，所述船舶水线信息包括首尾驻点空间坐标、水线拟合曲线函数；
17.步骤102:利用网格信息进行静水力计算具体包括，基于船舶物面网格，计算船舶静水力参数，检查静水力参数数值结果和物理船舶参数的误差，所述静水力参数包括排水体积、浮心、重心、漂心、惯性矩、横摇固有频率和船舶湿表面面积。
18.具体而言，所述步骤2具体包括：
19.步骤201：将浮体湿表面离散为若干四边形和三角形单元，在每一单元上，取单元节点坐标均值为中心，在面元中点对偶极强度作泰勒展开并保留一阶导数项，对源强作泰勒展开只保留一阶导数项；
20.步骤202:引入场点的切向一阶导数来封闭方程组，构成关于偶极强度、偶强的一阶导数为未知数，源强为已知变量的线性代数方程组；
21.步骤203:利用阻尼区来实现有航速船舶远方辐射问题，并且自由面条件采用积分格式；
22.步骤204:本方法求解的扰动速度势中包含辐射速度势和绕射速度势；
23.步骤205:扰动速度势求解的边界积分方程，采用矩阵求逆法进行时域步进求解。
24.具体而言，所述步骤3具体包括：
25.步骤301:通过伯努利方程获得船体表面压力，并在船体湿表面直接积分获得一阶波浪力，然后由6自由度运动方程获得船舶加速度，进而获得船舶的波频运动；
26.步骤302:在水平和首摇方向附加软弹簧，对船舶进行约束，防止数值发散。
27.具体而言，所述步骤4具体包括：
28.步骤401：利用步骤3获得的新时刻船舶波频运动更新扰动速度势边值条件，并再次使用泰勒展开边界元法求解扰动速度势，进而获得新时刻的船舶波频运动。重复步骤3和步骤4n
t
步，获得n
t
步带漂角船舶所受的二阶波浪漂移载荷。
29.具体而言，所述步骤5具体包括：
30.步骤501：将半个波浪周期内的二阶波漂载荷进行最小二乘拟合，并且直接带入mmg模型，获得下一时刻的船舶运动速度和空间位置信息。
31.本发明的有益效果在于：本发明提供了基于泰勒展开边界元法的波浪中船舶回转运动预报方法，该方法可以有效提高带漂角船舶二阶波浪漂移载荷计算精度和船舶回转轨迹数值模拟精度。
附图说明
32.图1是本发明船舶回转运动数值方法流程图；
33.图2是本发明船体和自由液面网格示意图；
34.图3是本发明空间固定坐标系和随船平动坐标系示意图；
35.图4是本发明kvlcc2船在迎浪条件下+35
°
回转数值结果示意图。
具体实施方式
36.为使本发明实施例的目的、技术方案和优点更加清楚，下面将对本发明实施例的技术方案进行清楚、完整地描述。
37.步骤1：读取网格文件，自动提取船舶水线信息和自动生成符合边界元方法的水面离散网格，所述船舶水线信息包括首尾驻点空间坐标、水线拟合曲线函数。基于船舶物面网格，计算船舶静水力参数，检查静水力参数数值结果和物理船舶参数的误差，所述静水力参数包括排水体积、浮心、重心、漂心、惯性矩、横摇固有频率和船舶湿表面面积。
38.步骤2：本发明利用泰勒展开边界元方法求解各速度势成分及其空间一阶、二阶导数。泰勒展开边界元方法核心思想是基于格林第三公式形成的边界积分方程进行数值离散求解的方法。对于三维问题，将浮体湿表面离散为若干四边形或三角形，在每一单元上，取单元节点坐标均值为中心，在面元中点对偶极强度作泰勒展开并保留一阶导数项，对源强作泰勒展开只保留一阶导数项。并引入场点的切向一阶导数来封闭方程组，从而构成了关于偶极强度、偶极的一阶导数为未知数，源强为已知变量的线性代数方程组。其中偶强切向一阶导数的影响系数包含主值。对于任意面元i可得到如下简化的一阶泰勒展开边界元方法的离散方程组，i＝1,2,
…
,n,
[0039][0040]
上式各矩阵中元素表达式：
[0041][0042][0043][0044]
式中：i和j表示面元编号。以矩阵中某一元素表达式为例做一解释，如：
[0045][0046]
该方法可同时求解速度势及两个相互正交的切向方向导数。在利用物面法向不可穿透条件即构成了当地局部坐标系的速度场。可实现速度场在局部坐标下及大地坐标系下的转换。数值结果表明，该方法可明显改善流域边界拐角处的切向诱导速度的计算精度。
[0047]
总速度势可被分解为基础速度势、定常速度势和非定常速度势，其中基础速度势可分解为来流速度势和叠模速度势，非定常速度势可分解为入射势和扰动势。即：
[0048][0049]
由于定常兴波对低航速船舶影响较小，因此，定常速度势可被忽略。非定常速度势定解问题的物面条件涉及叠模速度势的影响，即mj项。基于db假设，叠模速度势的边值定解问题为：
[0050][0051]
其中sh为船体平均湿表面，sf为静水面。在求得叠模速度势φ后，mj项表示为：
[0052][0053][0054]
其中其中为船体湿表面上各点位置矢径。
[0055]
本发明利用阻尼区处理带航速船舶远方辐射条件。如下为非定常扰动速度势边值问题:
[0056][0057]
其中线性化入射速度势已知，μ＝3μ0(l-l0)2/l
2-2μ0(l-l0)3/l3，
[0058]
l是距离坐标系原点的辐射距离，阻尼区起始边缘位于l＝l0，μ0和l分别为阻尼区阻尼强度及阻尼区长度。(n1,n2,n3)＝nj,ηj(j＝1-6)表示船舶的6自由度运动位移。
[0059]
步骤3：船舶所受的一阶波浪力包括回复力，f-k力，辐射力和绕射力。回复力通过回复力系数矩阵计算。在求解完扰动速度势后，利用伯努利方程在平均湿表面上积分，可获得一阶扰动波浪力，即f-k力，辐射力和绕射力。
[0060][0061]
在获得船舶所受的一阶波浪力后，根据牛顿第二定律可以求得船舶的波频运动：
[0062][0063]
当速度势及其空间一阶、二阶导数和船舶六自由度运动时历信号计算完成后，利
用近场压力积分公式，保留二阶压力载荷项，如下式所示为二阶波浪漂移力和力矩：
[0064][0065][0066]
其中表示船舶表面任一点的位移。表示船舶静止时的法向量。
[0067][0068][0069]
步骤5：本方法利用如下式对二阶波浪漂移力和力矩进行最小二乘拟合。
[0070]
y＝a0+a
1 cosωt+b
1 sinωt+a2cos2ωt+b
2 sin2ωt
[0071]
在图3随船平动坐标系下建立3自由度mmg操纵运动方程。计算船舶运动过程中，假定船舶为运动的刚体，其形状和内部质量分布不随运动而改变，并且船舶的纵向速度为u，横向速度为v，合速度为首摇角速度为r，船舶的基本操纵方程如下式：
[0072][0073][0074][0075]
式中m表示船舶的质量，i
zz
表示船舶绕z轴的惯性矩，xg为重心纵向坐标，方程右端表示船舶所受到的水动力和力矩。在操纵运动过程中，作用在船体上的水动力包括惯性力和非惯性力两部分。在本研究中，水动力和力矩可以用下式表示：
[0076][0077][0078]
[0079]
其中m
x
,my,j
zz
分别表示为船舶纵向、横向附加质量以及首摇方向的附加惯性矩。h,r,p分别表示作用于船体，舵和螺旋桨上的低频水动力。本模型忽略螺旋桨转动引起的横向力及其对重心的力矩。w表示作用于船体上的二阶波浪漂移载荷。因此，利用如下式的船舶操纵运动方程，求解下一时刻船舶的加速度，速度，艏向角和空间位置。
[0080][0081][0082][0083]
其中作用于船体上的低频水动力可被表示为：
[0084]
xh＝(1/2)ρl
pp
du2x'h(v',r')
[0085]
yh＝(1/2)ρl
pp
du2y'h(v',r')
[0086]
nh＝(1/2)ρl
2pp
du2n'h(v',r')
[0087]
v'和r'表示无因次化的横向速度和首摇角速度，而x'h,y'h和n'h可被表示为下式：
[0088]
x'h(v',r')＝-r'0+x'
vv
v'2+x'
vr
v'r'+x'
rr
r'2+x'
vvvv
v'4[0089]
y'h(v',r')＝y'vv'+y'rr'+y'
vvv
v'3+y'
vvr
v'2r'+y'
vrr
v'r'2+y'
rrr
r'3[0090]
n'h(v',r')＝n'vv'+n'rr'+n'
vvv
v'3+n'
vvr
v'2r'+n'
vrr
v'r'2+n'
rrr
r'3[0091]
作用于螺线桨上的低频水动力可被表示为：
[0092]
x
p
＝(1-t
p
)t
[0093]
其中t
p
是推力减额系数，t是螺旋桨推力，如下：
[0094]
t＝ρn
p2dp4kt
(j
p
)
[0095]
其中k
t
,d
p
和n
p
分别表示推力系数、螺旋桨直径和转速。k
t
可由螺旋桨进速系数j
p
表达为k
t
(j
p
)＝k2j
p2
+k1j
p
+k0，w
p
为螺旋桨的伴流系数，可被表达为：
[0096]
(1-w
p
)/(1-w
p0
)＝1+{1-exp(-c1|β
p
|)}(c
2-1)
[0097]
c1和c2为实验常数，β
p
＝β-x'
p
r'＝tan-1
(-v/u)-x'
p
r'。
[0098]
作用于舵上的低频水动力可被表示为：
[0099]
xr＝-(1-tr)fnsinδ
[0100]
yr＝-(1+αh)fncosδ
[0101]
nr＝-(xr+αhxh)fncosδ
[0102]
其中fn,δ和xr分别表示为舵的正压力、舵角和舵的纵向位置，αh,tr和xh分别表示为舵力修正因子、舵阻力减额系数和附加横向力的纵坐标。fn＝(1/2)ρaru
2rfα
sinαr，ar和f
α
[0103]
表示舵面积和舵的升力系数，ur和vr表示流入舵的纵向和横向速度。
[0104][0105]
vr＝uγrβr[0106]
βr＝β-l'
r r'
[0107]
其中ε和κ是实验常数，l'r和γr分别表示为舵无因次化纵向位置和整流系数。
[0108]
步骤6：通过步骤5获得的新时刻船舶运动信息，更新叠模速度势和扰动速度势的边值条件，进而计算下一时刻的船舶波频运动和二阶波浪漂移载荷。
[0109]
本发明提供了一种基于泰勒展开边界元法的波浪中船舶回转运动预报方法。该方法使用泰勒展开边界元法进行船舶耐波性计算，克服了传统常值面元法和高阶边界元法在求解非光滑边界处诱导速度精度低的缺点，有效提高了全浪向下带漂角船舶二阶波浪漂移力和力矩的计算精度，然后使用双时间尺度模型进行波浪中船舶回转运动计算。包括如下步骤：步骤1:读取网格文件，利用网格信息进行静水力计算。步骤2：通过时域泰勒展开边界元法计算船舶初始状态下的叠模速度势、扰动速度势及其空间一阶、二阶导数。步骤3：直接在船体湿表面压力积分获得一阶力，进而通过船舶6自由度运动方程求得船舶波频运动，并且通过近场积分方法直接获得船舶所受的二阶波浪漂移力和力矩。步骤4：利用步骤3获得的船舶6自由度波频运动更新扰动势边值问题，并且再次利用时域泰勒展开边界元法对扰动势进行求解，获得新时刻的船舶6自由度波频运动。步骤2到步骤4称为船舶耐波性问题，并且对耐波性问题进行nt步求解。步骤5：将获得的二阶波浪漂移力和力矩进行最小二乘拟合，并带入mmg操纵运动模型，计算下一时刻的船舶运动信息。利用本发明能够预报船舶在规则波中的运动rao，波漂载荷rao和回转轨迹。

技术特征：

1.一种基于泰勒展开边界元法的波浪中船舶回转运动预报方法，其特征在于，包括以下步骤：步骤1:读取网格文件，利用网格信息进行静水力计算；步骤2:通过时域泰勒展开边界元法计算船舶初始状态下的叠模速度势、扰动速度势及其空间一阶、二阶导数；步骤3:直接在船体湿表面压力积分获得一阶力，进而通过船舶6自由度运动方程求得船舶的波频运动，并且通过近场积分方法求得带漂角船舶所受的二阶波浪漂移力和力矩；步骤4：利用步骤3获得的船舶6自由度波频运动更新扰动势边值问题，并且再次利用时域泰勒展开边界元法对扰动势进行求解，获得新时刻的船舶6自由度波频运动，进而求得新时刻带漂角船舶所受的二阶波浪漂移力和力矩，步骤2到步骤4称为船舶耐波性问题，并且对耐波性问题进行n
t
步求解；步骤5：将获得的二阶波浪漂移力和力矩进行最小二乘拟合，并带入mmg操纵运动模型，计算下一时刻的船舶运动信息；步骤6：根据新的船舶运动信息重新计算叠模速度势和扰动速度势；步骤7:更新船舶的6自由度运动和带漂角船舶所受的二阶波浪漂移力和力矩，从而计算出下一时刻的船舶运动信息。2.根据权利要求1所述的基于泰勒展开边界元法的波浪中船舶回转运动预报方法，其特征在于，所述步骤1具体包括：步骤101:读取网格文件，自动提取船舶水线信息和自动生成符合边界元方法的水面离散网格，所述船舶水线信息包括首尾驻点空间坐标、水线拟合曲线函数；步骤102:利用网格信息进行静水力计算具体包括，基于船舶物面网格，计算船舶静水力参数，检查静水力参数数值结果和物理船舶参数的误差，所述静水力参数包括排水体积、浮心、重心、漂心、惯性矩、横摇固有频率和船舶湿表面面积。3.根据权利要求1所述的基于泰勒展开边界元法的波浪中船舶回转运动预报方法，其特征在于，所述步骤2具体包括：步骤201：将浮体湿表面离散为若干四边形和三角形单元，在每一单元上，取单元节点坐标均值为中心，在面元中点对偶极强度作泰勒展开并保留一阶导数项，对源强作泰勒展开只保留一阶导数项；步骤202:引入场点的切向一阶导数来封闭方程组，构成关于偶极强度、偶强的一阶导数为未知数，源强为已知变量的线性代数方程组；步骤203:利用阻尼区来实现有航速船舶远方辐射问题，并且自由面条件采用积分格式；步骤204:本方法求解的扰动速度势中包含辐射速度势和绕射速度势；步骤205:扰动速度势求解的边界积分方程，采用矩阵求逆法进行时域步进求解。4.根据权利要求1所述的基于泰勒展开边界元法的波浪中船舶回转运动预报方法，其特征在于，所述步骤3具体包括：步骤301:通过伯努利方程获得船体表面压力，并在船体湿表面直接积分获得一阶波浪力，然后由6自由度运动方程获得船舶加速度，进而获得船舶的波频运动；步骤302:在水平和首摇方向附加软弹簧，对船舶进行约束，防止数值发散。
5.根据权利要求1所述的基于泰勒展开边界元法的波浪中船舶回转运动预报方法，其特征在于，所述步骤4具体包括：步骤401：利用步骤3获得的新时刻船舶波频运动更新扰动速度势边值条件，并再次使用泰勒展开边界元法求解扰动速度势，进而获得新时刻的船舶波频运动，重复步骤3和步骤4n
t
步，获得n
t
步带漂角船舶所受的二阶波浪漂移载荷。6.根据权利要求1所述的基于泰勒展开边界元法的波浪中船舶回转运动预报方法，其特征在于，所述步骤5具体包括：步骤501：将半个波浪周期内的二阶波漂载荷进行最小二乘拟合，并且直接带入mmg模型，获得下一时刻的船舶运动速度和空间位置信息。

技术总结

本发明提供一种基于泰勒展开边界元法的波浪中船舶回转运动预报方法，步骤1：进行静水力计算；步骤2：通过时域泰勒展开边界元法计算船舶初始状态下的叠模速度势、扰动速度势及其空间一阶、二阶导数；步骤3：求得带漂角船舶所受的二阶波浪漂移力和力矩；步骤4：求得新时刻带漂角船舶所受的二阶波浪漂移力和力矩，并且对耐波性问题进行N