一种双浮体平台波浪激励力及运动响应预报方法

摘要

本发明属于多浮体海洋结构物的运动响应预报技术领域，具体涉及一种双浮体平台波浪激励力及运动响应预报方法。本发明构建了一种新的格林函数，该函数引入了波浪调节因子，增加了函数对波浪的敏感性，提高了算法的精度。本发明将新的格林函数应用于边界积分方程，通过对边界积分方程的离散处理，并与快速多极子技术结合，形成了一种新的快速多极子边界元方法，从而改进了势流问题中海洋结构物的应用范围，提高了计算效率。与传统的边界元方法相比，本发明对波浪较为敏感，且计算精度高、速度快，能够高效精确的计算多浮体海上结构物的水动力系数，在大型海洋结构物的水动力计算方面具有较大的优势。

权利要求

1.一种双浮体平台波浪激励力及运动响应预报方法，其特征在于，包括以下步骤：步骤1：输入待预报的双浮体平台及其运动状态j，获取双浮体平台的整体边界区域S0、第一浮体的边界区域SbA和第二浮体的边界区域SbB；输入环境参数，包括波浪频率ω和波高H；其中，j＝{1,2,3,4,5,6}分别表示横荡、纵荡、垂荡、横摇、纵摇、艏摇运动状态；

步骤2：计算双浮体平台处于第j运动状态下，满足关于速度势的边界积分方程的解pj、qj；

所述的关于速度势的边界积分方程为：

其中，pj、qj为双浮体平台的整体边界区域S0的两个点，点pj的坐标为(x1j,y1j,z1j)，点qj的坐标为(x2j,y2j,z2j)；α(qj)表示点qj的边界光滑度； 表示格林函数G(pj,qj)关于点pj垂线方向的偏导数；φ()为速度势； 表示点qj的速度势φ(qj)关于点qj垂线方向的偏导数；

格林函数G(pj,qj)的表达式为：

其中，fk为波浪调节因子；h、λ、μ均为常数；

J0()为贝塞尔函数；

速度势φ()的表达式为：

‑iωt ‑iωt ‑iωt

φ(x,y,z)＝Re[φI(x,y,z)e +φD(x,y,z)e +φR(x,y,z)e ]其中，φI(x,y,z)表示入射势，φD(x,y,z)表示绕射势，φR(x,y,z)表示辐射势；

步骤3：计算双浮体平台处于第j运动状态下的水动力系数；

其中，kAj为双浮体平台处于运动状态j时，第一浮体的运动方向；kBj为双浮体平台处于运动状态j时，第二浮体的运动方向；

步骤4：计算双浮体平台的第一浮体受到的波浪激励力的预报值fexA、第二浮体受到的波浪激励力的预报值fexB；

步骤5：计算双浮体平台的第一浮体运动响应的预报值RAOA、第二浮体运动响应的预报值RAOB；

ξA和ξB根据运动方程求解得到：

2 2

(‑ω(MAj+aAA)‑iωbAA+CA)ξA+(‑ωaAB‑iωbAB)ξB＝fexA+FLAj

2 2

(‑ωaBA‑iωbBA)ξA+(‑ω(MBj+aBB)‑iωbBB+CB)ξB＝fexB+FLBj其中，MAj为双浮体平台处于运动状态j时，第一浮体的物体质量；MBj为双浮体平台处于运动状态j时，第二浮体的物体质量；CA为第一浮体的静水回复力；CB为第二浮体的静水回复力；当j＝{1,2,3}，即双浮体平台处于横荡、纵荡或垂荡运动状态时，FLAj表示第一浮体与第二浮体之间的连接装置对第一浮体的作用力，FLBj表示第一浮体与第二浮体之间的连接装置对第二浮体的作用力；当j＝{4,5,6}，即双浮体平台处于横摇、纵摇或艏摇运动状态时，FLAj表示第一浮体与第二浮体之间的连接装置对第一浮体的作用力矩，FLBj表示第一浮体与第二浮体之间的连接装置对第二浮体的作用力矩；FLAj与FLBj大小相等，方向相反。

说明书

一种双浮体平台波浪激励力及运动响应预报方法

技术领域

[0001] 本发明属于多浮体海洋结构物的运动响应预报技术领域，具体涉及一种双浮体平台波浪激励力及运动响应预报方法。

背景技术

[0002] 在现有的数值方法中，边界元法(BEM)相对于其他数值算法有显著的优点，在选定格林函数后，边界元法将对整个域的求解转移到边界上，从而只需要在海洋结构的表面进行网格划分和未知量的设定，减少了计算存储量及计算时间，前期的准备也相对较少。传统的边界元法具有很多优点，但对边界积分方程(BIE)进行离散后，产生线性方程组的系数矩阵导致方程组的求解需要大量的计算时间与计算量，且均以未知数的平方量级增加。当未知量比较少时，BEM的计算效率较高。但针对大型海洋结构物，边界元的个数非常大，计算量和存储空间倍增，对计算机的能力需求也随之增加，使得边界元法优势不能充分利用。

[0003] 而应用边界元法对多浮体海洋结构物进行水动力计算时也需要更多的未知单元，为克服传统边界元法针对大型海洋结构进行水动力计算的缺点，将快速多极子技术应用于边界元方法，即快速多极子边界元法，以解决多浮体海洋结构物水动力问题，该方法在提高效率的基础上又降低计算存储量。

[0004] 在使用边界元方对海洋结构的水动力性能进行分析时，随着结构物尺寸的加大，未知量的个数也逐渐增加，导致该方法在进行数值计算时需要的内存及计算量大幅增加，从而边界元方法针对结构物进行数值计算的优势被抵消。快速多极子技术可解决大量粒子间相互作用的计算问题，将其与边界元法的有效结合解决了关于大型结构物的工程问题，充分提高了计算效率。

[0005] 经检索发现，公布号为CN109344531 A的发明专利文件中公开了预报多浮体结构物波漂载荷的三维频域数值方法，其特征是：读取网格文件，利用网格信息进行船舶静水力计算；计算简单格林函数的影响系数矩阵；计算复杂频域格林函数的影响系数矩阵；应用泰勒展开边界元方法求解各单位辐射速度势；多浮体结构的各单体水动力系数求解；多浮体运动方程求解；多浮体结构的整体、各单体的波浪力和波浪载荷求解。本发明所涉及的一种基于改进满足自由面边界条件的格林函数的多浮体海洋结构物数值计算方法与公布号为CN109344531 A的发明专利文件中公开的预报多浮体结构数值计算方法存在以下不同：

[0006] 1、格林函数的不同。公布号为CN109344531A的发明专利文件中公开的格林函数分别为简单格林函数与复杂频域格林函数。本发明采用的格林函数不仅满足自由面边界条件，而且增加了波浪调节因子。采用本发明所述的格林函数，增加了算法对波浪的敏感性，提高了算法的精度。

[0007] 2、边界元的计算方式不同。公布号为CN109344531A的发明专利文件中公开的多浮体水动力系数求解方法应用泰勒展开边界元法求解。本发明采用的方法为快速多极子技术与边界元结合的快速多极子边界元法。采用本发明所述结合改进自由面格林函数的快速多极子边界元法，提高了计算的效率，便于针对大型海洋结构物及多浮体海洋结构物的计算。

发明内容

[0008] 本发明的目的在于提供一种双浮体平台波浪激励力及运动响应预报方法。

[0009] 本发明的目的通过如下技术方案来实现：包括以下步骤：

[0010] 步骤1：输入待预报的双浮体平台及其运动状态j，获取双浮体平台的整体边界区域S0、第一浮体的边界区域SbA和第二浮体的边界区域SbB；输入环境参数，包括波浪频率ω和波高H；其中，j＝{1,2,3,4,5,6}分别表示横荡、纵荡、垂荡、横摇、纵摇、艏摇运动状态；

[0011] 步骤2：计算双浮体平台处于第j运动状态下，满足关于速度势的边界积分方程的解pj、qj；

[0012] 所述的关于速度势的边界积分方程为：

[0013]

[0014] 其中，pj、qj为双浮体平台的整体边界区域S0的两个点，点pj的坐标为(x1j,y1j,z1j)，点qj的坐标为(x2j,y2j,z2j)；α(qj)表示点qj的边界光滑度； 表示格林函数G(pj,qj)关于点pj垂线方向的偏导数；φ()为速度势； 表示点qj的速度势φ(qj)关于点qj垂线方向的偏导数；

[0015] 格林函数G(pj,qj)的表达式为：

[0016]

[0017] 其中，fk为波浪调节因子；h、λ、μ均为常数；J0()为贝塞尔函数；

[0018] 速度势φ()的表达式为：

[0019] φ(x,y,z)＝Re[φI(x,y,z)e‑iωt+φD(x,y,z)e‑iωt+φR(x,y,z)e‑iωt][0020] 其中，φI(x,y,z)表示入射势，φD(x,y,z)表示绕射势，φR(x,y,z)表示辐射势；

[0021] 步骤3：计算双浮体平台处于第j运动状态下的水动力系数；

[0022]

[0023]

[0024]

[0025]

[0026] 其中，kAj为双浮体平台处于运动状态j时，第一浮体的运动方向；kBj为双浮体平台处于运动状态j时，第二浮体的运动方向；

[0027] 步骤4：计算双浮体平台的第一浮体受到的波浪激励力的预报值fexA、第二浮体受到的波浪激励力的预报值fexB；

[0028]

[0029]

[0030] 步骤5：计算双浮体平台的第一浮体运动响应的预报值RAOA、第二浮体运动响应的预报值RAOB；

[0031]

[0032]

[0033] ξA和ξB根据运动方程求解得到：

[0034] (‑ω2(MAj+aAA)‑iωbAA+CA)ξA+(‑ω2aAB‑iωbAB)ξB＝fexA+FLAj[0035] (‑ω2aBA‑iωbBA)ξA+(‑ω2(MBj+aBB)‑iωbBB+CB)ξB＝fexB+FLBj[0036] 其中，MAj为双浮体平台处于运动状态j时，第一浮体的物体质量；MBj为双浮体平台处于运动状态j时，第二浮体的物体质量；CA为第一浮体的静水回复力；CB为第二浮体的静水回复力；当j＝{1,2,3}，即双浮体平台处于横荡、纵荡或垂荡运动状态时，FLAj表示第一浮体与第二浮体之间的连接装置对第一浮体的作用力，FLBj表示第一浮体与第二浮体之间的连接装置对第二浮体的作用力；当j＝{4,5,6}，即双浮体平台处于横摇、纵摇或艏摇运动状态时，FLAj表示第一浮体与第二浮体之间的连接装置对第一浮体的作用力矩，FLBj表示第一浮体与第二浮体之间的连接装置对第二浮体的作用力矩；FLAj与FLBj大小相等，方向相反。

[0037] 本发明的有益效果在于：

[0038] 本发明构建了一种新的格林函数，该函数引入了波浪调节因子，增加了函数对波浪的敏感性，提高了算法的精度。本发明将新的格林函数应用于边界积分方程，通过对边界积分方程的离散处理，并与快速多极子技术结合，形成了一种新的快速多极子边界元方法，从而改进了势流问题中海洋结构物的应用范围，提高了计算效率。与传统的边界元方法相比，本发明对波浪较为敏感，且计算精度高、速度快，能够高效精确的计算多浮体海上结构物的水动力系数，在大型海洋结构物的水动力计算方面具有较大的优势。

附图说明

[0039] 图1为快速多极子边界元法与边界元法存储量的变化图。

[0040] 图2为快速多极子边界元法与边界元法计算时间的变化图。

[0041] 图3为辐射阻尼对比图。

[0042] 图4为附加质量对比图。

[0043] 图5为几何构型对波浪激励力的影响图。

[0044] 图6为几何构型对俘获宽度比的影响图。

具体实施方式

[0045] 下面结合附图对本发明做进一步描述。

[0046] 本发明涉及多浮体海洋结构物及大型海洋结构物数值计算，属于势流理论的计算方法分析领域。

[0047] 本发明提供了一种双浮体平台波浪激励力及运动响应预报方法，解决了具有波敏性的三维频域水动力性能计算的数值问题。本发明对波浪较为敏感，且计算精度高、速度快，与传统的边界元方法相比，在大型海洋结构物的水动力计算方面具有较大的优势。

[0048] 一种双浮体平台波浪激励力及运动响应预报方法，包括以下步骤：

[0049] 步骤1：输入待预报的双浮体平台及其运动状态j，获取双浮体平台的整体边界区域S0、第一浮体的边界区域SbA和第二浮体的边界区域SbB；输入环境参数，包括波浪频率ω和波高H；其中，j＝{1,2,3,4,5,6}分别表示横荡、纵荡、垂荡、横摇、纵摇、艏摇运动状态；

[0050] 步骤2：计算双浮体平台处于第j运动状态下，满足关于速度势的边界积分方程的解pj、qj；

[0051] 所述的关于速度势的边界积分方程为：

[0052]

[0053] 其中，pj、qj为双浮体平台的整体边界区域S0的两个点，点pj的坐标为(x1j,y1j,z1j)，点qj的坐标为(x2j,y2j,z2j)；α(qj)表示点qj的边界光滑度； 表示格林函数G(pj,qj)关于点pj垂线方向的偏导数；φ()为速度势； 表示点qj的速度势φ(qj)关于点qj垂线方向的偏导数；

[0054] 格林函数G(pj,qj)的表达式为：

[0055]

[0056] 其中，fk为波浪调节因子；h、λ、μ均为常数；J0()为贝塞尔函数；

[0057] 速度势φ()的表达式为：

[0058] φ(x,y,z)＝Re[φI(x,y,z)e‑iωt+φD(x,y,z)e‑iωt+φR(x,y,z)e‑iωt][0059] 其中，φI(x,y,z)表示入射势，φD(x,y,z)表示绕射势，φR(x,y,z)表示辐射势；

[0060] 步骤3：计算双浮体平台处于第j运动状态下的水动力系数；

[0061]

[0062]

[0063]

[0064]

[0065] 其中，kAj为双浮体平台处于运动状态j时，第一浮体的运动方向；kBj为双浮体平台处于运动状态j时，第二浮体的运动方向；

[0066] 步骤4：计算双浮体平台的第一浮体受到的波浪激励力的预报值fexA、第二浮体受到的波浪激励力的预报值fexB；

[0067]

[0068]

[0069] 步骤5：计算双浮体平台的第一浮体运动响应的预报值RAOA、第二浮体运动响应的预报值RAOB；

[0070]

[0071]

[0072] ξA和ξB根据运动方程求解得到：

[0073] (‑ω2(MAj+aAA)‑iωbAA+CA)ξA+(‑ω2aAB‑iωbAB)ξB＝fexA+FLAj[0074] (‑ω2aBA‑iωbBA)ξA+(‑ω2(MBj+aBB)‑iωbBB+CB)ξB＝fexB+FLBj[0075] 其中，MAj为双浮体平台处于运动状态j时，第一浮体的物体质量；MBj为双浮体平台处于运动状态j时，第二浮体的物体质量；CA为第一浮体的静水回复力；CB为第二浮体的静水回复力；当j＝{1,2,3}，即双浮体平台处于横荡、纵荡或垂荡运动状态时，FLAj表示第一浮体与第二浮体之间的连接装置对第一浮体的作用力，FLBj表示第一浮体与第二浮体之间的连接装置对第二浮体的作用力；当j＝{4,5,6}，即双浮体平台处于横摇、纵摇或艏摇运动状态时，FLAj表示第一浮体与第二浮体之间的连接装置对第一浮体的作用力矩，FLBj表示第一浮体与第二浮体之间的连接装置对第二浮体的作用力矩；FLAj与FLBj大小相等，方向相反。

[0076] 本发明所使用的格林函数不仅满足自由面边界条件，而且添加了以波数为变量的波浪调节因子fk，增加了对波浪的敏感度。本发明所使用的一种新的格林函数的表达式为：

[0077]

[0078] 其中波数线性函数fk，可针对不同周期下的波浪特性，对波数线性函数中波浪调节因子进行调节，使得其满足不同的波浪特征。

[0079] 本发明具体包括以下五个步骤：

[0080] [1]读取多浮体边界离散后的网格文件，利用多极子技术进行边界积分方程的计算；

[0081] [2]将离散后的边界积分中包含以波数为变量的波浪调节因子的格林函数展开成源点和场点的乘积；

[0082] [3]对含有波数线性函数的格林函数通过展开后的调和函数进行多极矩及多极展开的计算；

[0083] [4]根据离散后源点及场点的分布将其划分为近场及远场两部分，针对近场采用传统边界元法直接计算，而远场则采用新的快速多极子技术进行计算，最后将近场部分与远场部分相加，从而求得速度势；

[0084] [5]依据速度势计算结果得到多浮体水动力系数，分析水动力性能。

[0085] 本发明建了一种新的格林函数，该函数引入了波浪调节因子，增加了函数对波浪的敏感性，提高了算法的精度。本发明将新的格林函数应用于边界积分方程，通过对边界积分方程的离散处理，并与快速多极子技术结合，形成了一种新的的快速多极子边界元方法。从而改进了势流问题中海洋结构物的应用范围，提高了计算效率。利用本发明提出的方法能够高效精确的计算多浮体海上结构物的水动力系数，为装置的设计提供了准确的计算依据。

[0086] 假设流体为无粘、无旋、不可压缩的理想流体，流场内流体速度势Φ满足Laplace方程 对于任意空间Ω中两个二阶连续可微函数φ和G，使得格林公式成立 格林公式需满足的条件有：自由
面条件 远方条件 底部条件 因此关于
速度势的边界积分方程可表示为 其
‑iωt ‑iωt ‑iωt
中速度势为Φ(x,y,z)＝Re(φ(x,y,z)e )＝Re[φI(x,y,z)e +φD(x,y,z)e +φR‑iωt
(x,y,z)e ]

[0087] 其中φI(x,y,z)表示入射势，φD(x,y,z)表示绕射势，φR(x,y,z)表示辐射势，根据控制方程及不同的物面条件可分别求出速度势中的绕射势和辐射势。

[0088] 针对多体边界方程通过引入八个点的插值函数，并将系数应用高斯积分进行计算，将整体坐标系换为局部坐标系，则通过二次型离散化处理可得离散后的边界积分方程为

[0089]

[0090]

[0091]

[0092]

[0093] …

[0094] 式中ui*、ui*′、ui*″、ui*″′表示格林函数， 为格林函数的法向偏导。

[0095] 求解离散后的边界积分方程采用边界元法与快速多极子技术相结合进行求解，其中远场部分采用上述的快速多极子技术计算，近场部分采用边界元方法进行计算。首先将离散后的点进行编号及层级划分，形成树结构；其次对含有以波数为变量调节因子的格林函数展开成源点与场点的乘积形式；最后采用M2M、M2L、L2L三种多极子展开方法对改进满足自由面的格林函数进行局部展开系数及多极动量矩的计算。

[0096] ∫SG(x,y)q(x)dS(x)＝∫S_nearG(x,y)q(x)dS(x)+∫S_farG(x,y)q(x)dS(x)[0097] 根据本文涉及发明的速度势计算方法即可求得多浮体运动时自身产生的水动力系数即其他浮体对自身影响的水动力系数。

[0098]

[0099]

[0100]

[0101]

[0102] aAA、bAA、aBB及bBB表示为浮体在运动时产生水动力系数，aAB、bAB、aBA及bBA表示在其他浮体作用下产生影响的水动力系数。

[0103] 则波浪激励力及运动响应表示为如下四个表达式。

[0104] 波浪激励力：

[0105]

[0106]

[0107] 运动响应：

[0108] ξA可根据运动方程求得：

[0109] (‑ω2(MAj+aAA)‑iωbAA+CA)ξA+(‑ω2aAB‑iωbAB)ξB＝fexA+FLAj[0110] (‑ω2aBA‑iωbBA)ξA+(‑ω2(MBj+aBB)‑iωbBB+CB)ξB＝fexB+FLBj[0111] 其中，MAj为双浮体平台处于运动状态j时，第一浮体的物体质量；MBj为双浮体平台处于运动状态j时，第二浮体的物体质量；CA为第一浮体的静水回复力；CB为第二浮体的静水回复力；当j＝{1,2,3}，即双浮体平台处于横荡、纵荡或垂荡运动状态时，FLAj表示第一浮体与第二浮体之间的连接装置对第一浮体的作用力，FLBj表示第一浮体与第二浮体之间的连接装置对第二浮体的作用力；当j＝{4,5,6}，即双浮体平台处于横摇、纵摇或艏摇运动状态时，FLAj表示第一浮体与第二浮体之间的连接装置对第一浮体的作用力矩，FLBj表示第一浮体与第二浮体之间的连接装置对第二浮体的作用力矩；FLAj与FLBj大小相等，方向相反，不同的连接形式的FL表达式不同，对于铰连接或滑槽连接，如果不能限制物体某一方向的运动，则该项为零。

[0112] 如图1所示，本发明专利的数值计算方法与边界元法计算存储量对比可知该方法计算量小的特点。看出随着未知量个数的增加，两种方法的存储量均呈上升趋势，且明显看出传统边界元法的上升速度较快，需要的存储量较大，且未知量个数越来越大时，本发明专利的数值计算方法的存储量要小于传统边界元计算时所需的存储量，充分验证了本发明专利的数值计算方法可以减小计算存储量的优势。所述计算具有波敏性的三维频域水动力数值方法降低了对计算能力的需求。

[0113] 如图2所示，随着未知量逐渐增大，两种方法的计算时间均呈递增趋势。当未知量个数小于4000时，边界元法的计算时间较短，效率较高，而当未知量个数高于4000，传统边界元的计算时间要远远高于FMM‑BEM的计算时间，因此FMM‑BEM更适合大规模多浮体未知数多的物理模型。

[0114] 如图3及图4所示，本发明专利提出的一种具有波敏性的三维频域水动力数值方法的计算结果与传统边界元法的数值计算结果进行对比可知，两种方法计算附加质量及辐射阻尼的计算结果误差较小，充分验证了该数值计算方法精确性。

[0115] 如图5及图6所示，应用本发明专利的高效性及精确性对不同几何构型浮子的水动力性能进行分析，通过本发明专利计算可得锥底柱形浮子的峰值和宽度带宽均大于另外两种浮子，因此FMM‑BEM为水动力数值分析提供了有效的方法。

[0116] 以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。