CN109726423B - 圆柱阵列波浪力幅值波动特性、步长和包络线的获取方法 - Google Patents

Abstract

本发明提供一种圆柱阵列波浪力幅值波动特性、步长和包络线的获取方法，对排成一条直线的有限个坐底圆柱所组成的圆柱阵列建立坐标系并确定相关参数，基于相长相消理论分析任意圆柱及其上下游圆柱绕射波的干涉效果，得到波浪入射角等于零和不等于零时波动间距的描述模型、波动间距和峰谷点横坐标的表达式，确定波浪力幅值曲线计算步长，以及峰谷点的纵坐标；将各峰/谷点连线得到波浪力曲线在非捕获区中的上/下包络线。本发明还提供了基于绕射理论建立波浪力幅值曲线波动间距描述模型的方法。本发明证明了波动间距只与圆柱阵列中圆柱总数、标识圆柱位置的柱子编号和波浪入射角有关，可在保证精度的前提下减少工作量，缩短设计和评估周期。

Description

圆柱阵列波浪力幅值波动特性、步长和包络线的获取方法

技术领域

本发明涉及海洋工程领域，特别涉及由大数量穿透水面圆柱组成的圆柱阵列在波浪作用下，基于相长相消理论建立任意圆柱所受波浪力的波浪力幅值在非捕获区波动间距的描述模型、波浪力幅值曲线计算步长的确定方法以及在非捕获区波浪力幅值曲线包络线的确定方法；还提供基于绕射理论建立任意圆柱所受波浪力的波浪力幅值在非捕获区波动间距的描述模型的方法。

背景技术

占地球表面积71％的海洋中蕴藏着丰富的石油、天然气等可采资源和风能、波浪能等可供人类长期使用的可再生能源。随着经济发展对能源和资源需求的不断增加，在海洋中拓展生存空间并寻求各种物质和能源供给已经成为明显趋势。

无论是进行海洋资源开采，海上空间开发，还是实际利用海洋可再生能源，都需要发展作为载体的海上结构物。这其中有很重要的一类结构物，尽管其上部建筑有所不同，但它们的浮体/支撑结构均是由多个穿透水面的圆柱(即，圆柱阵列)组成。例如，海洋石油平台，跨海大桥，超大浮体，波浪发电阵列等。随着经济社会发展对海洋开发需求的不断扩大，海洋结构物整体尺度变得越来越大，作为上述海洋结构物浮体/支撑结构的圆柱阵列的规模也随之增大。圆柱阵列中圆柱的数目从最初的个位数增加到十几个、几十个、数百个，甚至可达到上千个。单排坐底圆柱阵列是圆柱阵列的典型型式之一。此处单排坐底圆柱阵列指的是：水平面与各圆柱相交所得圆形横截面的圆心在一条直线上、从水底不间断地直通水面并穿透水面延伸向上的圆柱组成的圆柱阵列。水中圆柱阵列所受到的波浪力是决定圆柱阵列设计方案、保证结构安全性的一个关键要素，为此，需要掌握波浪力幅值随无量纲波数变化的规律。

如图1所示的单排坐底圆柱阵列中任意圆柱上所受波浪力的幅值是随无量纲波数上下起伏、波动变化的。一般来说，单排较大数量(例如，数量大于9)有限圆柱阵列中单个柱子所受到波浪力的幅值随波数变化具有三个明显的特点：1)波浪力幅值随无量纲波数变化形成的波浪力曲线具有几个高耸的尖峰，这几个高耸尖峰所在区域本发明称为区域I(RegionI)；2)在区域I附近，该曲线有若干个逐渐变小的次级峰和谷，这些次级峰的高度均低于区域I的尖峰且波浪力曲线波动间距随无量纲波数的改变而发生变化，该区域本发明称为区域II(Region II)；3)在上述两个区域之外，很多地方具有非常规则的波动性，该区域本发明称为区域III(Region III)。上述三个区域示意图如图2所示。

区域I和区域II与near-trapping有关，本发明称这两个区域为“捕获有关区域”。而区域III，本发明称其为“非捕获区”。与捕获有关的区域I和区域II，国际上已有不少near-trapping的研究发表，对其理解较为深刻。对于非捕获区，目前还没有对其波动规律有深入研究，更缺少一种描述非捕获区波动间距的描述模型。本发明所述波动间距是指：波浪力幅值随无量纲波数变化形成的波浪力曲线上，相邻两个极大值点(或极小值点)横坐标之间的距离。在本发明中，也用“峰”或“谷”来指代极大值点或极小值点。对非捕获区波动间距进行描述可以提高设计水平从而有助于以较小的代价延长结构的疲劳寿命。这是因为，考察大量计算结果后发现，波浪力曲线区域III中相邻峰和谷数值的相对差有时较大，仅从有限数量的计算结果中就可发现，在区域III中这种相对差最大可达20％左右。

因此在实际计算波浪力的过程中，如果横坐标步长取得不够小，那么在区域III波浪力计算结果的误差就有可能达到20％甚至更大。对于极限载荷引起的“一次性”强度破坏问题，这可能影响不大，因为捕获区尖峰处的波浪力幅值远高于非捕获区域，非捕获区域相对较小的波浪力幅值的20％左右误差不会对结构的“一次性破坏”产生影响。但是，对于循环载荷引起的疲劳破坏来说，由于疲劳寿命的计算需要计入一定频率范围内的波浪力的综合贡献(而不仅仅如强度分析那样仅需考虑捕获区near-trapping频率下对应的最大值)，因此上述波浪力计算误差可能产生明显的负面影响。这是因为，在分析线性时不变系统疲劳寿命时，交变应力响应的谱密度函数等于输入的海浪谱密度乘以系统传递函数模的平方。常规海洋结构物弹性模态固有频率远高于波浪频率，因此将如图2所示的波浪力幅值的传递函数乘以某个系数即可以得到交变应力幅值的传递函数。如果波浪力传递函数在计算时由于步长选择不恰当产生较大误差，则交变应力幅值传递函数同样也会产生较大误差，那么平方后这个误差就会变得更大(例如，若传递函数的模误差10％，平方后误差就会增加到20％，若传递函数的模误差20％，平方后误差就会增加到36％)。从而会得到不准确的交变应力响应结果，进而影响疲劳寿命评估的准确性。考虑到，一般情况下圆柱阵列在设计时会使系统的near-trapping频率避开海浪能量较大的频段，这会使非捕获区的交变应力在对疲劳损伤的贡献中占据很大份额。因而，如图2所示的非捕获区波浪力的准确计算对于疲劳寿命的准确评估有重要意义。

综上所述，想要高效准确地获取非捕获区波浪力的前提是，掌握非捕获区波浪力曲线波动特性，而且获得能够事先准确预测非捕获区波浪力曲线波动间距的描述模型。对于占波浪力曲线大部分区域的非捕获区，即区域III(该区域对评估结构的疲劳寿命有实际意义)，目前对其波动特性还缺少深刻的认识，尚没有波浪力幅值随无量纲波数变化形成的波浪力曲线在非捕获区波动间距的描述模型来作为进行高效评估和设计的基础。

此外，目前尽管对区域I和区域II的研究较多，但如前所述，由于对区域III波动特性还缺乏深刻认识，所以为得到准确的波浪力曲线，通常采用改变计算步长反复试算的方法进行。计算步长如何选取，在计算开始之前无法定量估算，基本是通过猜测和尝试来确定初始计算步长以及进而修正计算步长。这一过程即使对非常有经验的专家来说，也是很繁琐和费时费力的。对于经验不足或完全没有经验的人，这一过程则是非常繁重、代价高昂的。

另外一个问题是，目前，在基于波浪力曲线来进行有关实际工程结构设计和评估时，为避免漏掉波浪力曲线的峰谷点，通常需要对波浪力曲线进行繁复大量的计算，这会引起时间成本过大和计算代价高昂等问题。而在方案评估、选型或初步设计时，这种做法会导致效率低下，大幅增加成本、虚耗时间，非常不经济。此时，如果能够准确快速掌握波浪力曲线波动的上下界限，就可以在避免繁复计算的情况下给出合理的初步解决方案，因此迫切需要一种能够快速准确获得波浪力曲线峰点和谷点的方法，从而解决上述问题。

发明内容

本发明的目的是要提供一种由大数量穿透水面圆柱组成的圆柱阵列在波浪作用下，基于相长相消理论建立任意圆柱所受波浪力的波浪力幅值在非捕获区波动间距的描述模型、波浪力幅值曲线计算步长的确定方法以及在非捕获区波浪力幅值曲线包络线的确定方法；还提供基于绕射理论建立任意圆柱所受波浪力的波浪力幅值在非捕获区波动间距的描述模型的方法。

特别地，本发明提供圆柱阵列波浪力幅值波动特性、步长和包络线的获取方法，包括如下步骤：

步骤100，将波浪力幅值随无量纲波数变化形成的波浪力曲线中的多个高耸尖峰所在区域作为区域I，将高耸尖峰附近低于所述高耸尖峰且曲线波动间距随无量纲波数改变而变化的次级峰和谷所在区域作为区域II，将排除区域I和区域II后的波浪力曲线作为区域III；

步骤200，建立由排成一条直线的有限个相同坐底圆柱所组成的圆柱阵列坐标系并确定相关参数，将无量纲波数转化为相邻两圆柱间距离和波长的比值，则各圆柱的绕射波发生相长干涉时波浪力曲线出现峰点，发生相消干涉时波浪力曲线出现谷点；对于波浪入射角等于零的情形，分析入射波传到任意某个圆柱发生绕射及入射波传到该圆柱下游最后一个圆柱发生绕射这两个路径的波程差，求出波浪力曲线区域III中任意峰点横坐标的初步表达式和相邻峰点横坐标之差的初步表达式，以及任意谷点横坐标的初步表达式和相邻谷点横坐标之差的初步表达式，进而可知相邻两个峰点或相邻两个谷点的横坐标之差相等，最后得到任意圆柱波浪力曲线在区域III中的初步波动间距表达式；

步骤300，在任意某个圆柱上游处，将该圆柱产生的向上游传播的左传绕射波和该圆柱下游各个圆柱产生的左传绕射波求和，利用汉克尔函数取渐近表达式进行简化，得到对波动间距问题来说与当前圆柱阵列等价、且位置与现有圆柱位置均不相同的两个位于该圆柱下游的等价圆柱的具体位置，将离该圆柱更远的那个等价圆柱的位置代入初步波动间距表达式、任意峰点横坐标的初步表达式和任意谷点横坐标的初步表达式进行修正，即可得到波浪入射角等于零时波浪力曲线区域III中波动间距的最终表达式、任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式；

步骤400，对于波浪入射角不等于零的情形，利用前述同样方法，首先分析入射波传到任意某个圆柱发生绕射及入射波传到该圆柱上游的阵列端部第一个圆柱发生绕射这两个路径的波程差，然后分析入射波传到任意某个圆柱发生绕射及入射波传到该圆柱下游的最后一个圆柱发生绕射这两个路径的波程差，即可分别得到任意圆柱波浪力曲线在区域III中由于该圆柱上游圆柱绕射波的作用导致的任意峰点横坐标的初步表达式一、任意谷点横坐标的峰点初步表达式一和初步波动间距表达式一，以及由于该圆柱下游圆柱绕射波的作用导致的任意峰点横坐标的初步表达式二、任意谷点横坐标的初步表达式二和初步波动间距表达式二；

步骤500，首先在任意某个圆柱的下游处，对该圆柱产生的向下游传播的右传绕射波和该圆柱上游各个圆柱产生的右传绕射波求和，利用前述同样的方法修正，得到波浪入射角不等于零时，任意峰点横坐标的最终表达式一、任意谷点横坐标的最终表达式一和最终波动间距表达式一；然后，在任意某个圆柱的上游处，对该圆柱产生的向上游传播的左传绕射波和该圆柱下游各个圆柱产生的左传绕射波求和，同样方法修正后得到任意峰点横坐标的最终表达式二、任意谷点横坐标的最终表达式二和最终波动间距表达式二；

步骤600，采用最终波动间距表达式一和最终波动间距表达式二给出的波动间距的更小者，就得到波浪入射角不等于零时任意圆柱波浪力曲线在区域III中最小波动间距的最终表达式；综合波浪入射角等于零时和不等于零时的波动间距的最终表达式，即可得到任意圆柱波浪力曲线在区域III中波动间距描述模型的最终表达式；

步骤700，根据入射角等于零和不等于零时波动间距的最终表达式计算出区域III中的最小波动间距，以该最小波动间距作为波浪力曲线在区域III中计算步长的上限，根据不同精度要求将最小波动间距除以相应的自然数作为区域III中计算步长的下限，从而得到区域III的计算步长；

根据所需精度要求将区域III的最小波动间距除以相应精度要求对应的自然数，得到区域II中的计算步长；根据所需精度要求将区域III的最小波动间距除以相应精度要求对应的自然数，得到区域I中的计算步长；

步骤800，采用最终表达式一和最终表达式二中横坐标的更小者作为波浪入射角不等于零时任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式，综合波浪入射角等于零时和不等于零时的最终表达式，即可得到任意圆柱波浪力曲线在区域III中任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式；根据横坐标的最终表达式得到任意峰点和谷点处的横坐标，从而可得对应的波数，求解线性方程组以获得速度势表达式中的未知绕射系数，进而可得到任意圆柱所受波浪力，将该波浪力无量纲化并取模，就得到波浪力曲线在区域III中任意峰点和谷点处的纵坐标；

根据上述获取的任意峰点和谷点处的横坐标以及任意峰点和谷点处的纵坐标，确定出每一个峰点和谷点的位置，将这些峰点连线，就得到波浪力曲线在区域III中的上包络线；将这些谷点连线，就得到波浪力曲线在区域III中的下包络线。

在本发明的一个实施方式中，在所述圆柱阵列坐标系中，平面入射波传播方向与圆柱阵列整体坐标系中x轴正方向的夹角称为波浪入射角，且整体坐标系的建立使波浪入射角小于等于90度；k为圆柱阵列中的任意一根圆柱的编号，编号k增大方向与圆柱阵列整体坐标系中x轴正方向保持一致。

在本发明的一个实施方式中，所述步骤200中的初步波动间距表达式、任意峰点横坐标的初步表达式和任意谷点横坐标的初步表达式的获取过程如下：

无量纲波数Kd/π可以改写为：Kd/π＝R/λ；

波浪力曲线中的第一个峰点R_p(1)和谷点R_v(1)对应的柱间距用下式表示：

其中，K是波数，R＝2d是相邻圆柱的柱间距，λ是波长，脚标中，p代表峰点，v代表谷点，(1)代表第一个峰或者谷；

令在区域III中任意第s和第s+1个峰点对应的柱间距是R_p(s)和R_p(s+1)，R_p(s+1)＝R_p(s)+δR_p，峰点出现的条件是发生相长干涉，则相邻的第s和第s+1个峰点对应的波程差应该分别是s倍和s+1倍的波长，由此可得相邻峰点横坐标之差的初步表达式为：

谷点出现的条件是发生相消干涉，对于第s和第s+1个谷值点对应的柱间距R_v(s)和R_v(s+1)＝R_v(s)+δR_v，相应的波程差分别等于2s-1倍和2s+1倍半波长，由此可得相邻谷点横坐标之差的初步表达式为：

由上述推导可知，相邻峰点和相邻谷点的横坐标之差相等，则初步波动间距表达式为：

N为圆柱阵列中的圆柱总数；此时波浪力曲线区域III中任意峰点横坐标的初步表达式和任意谷点横坐标的初步表达式为：

在本发明的一个实施方式中，获取波浪入射角等于零时波浪力曲线区域III中波动间距的最终表达式、任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式的过程如下：

在任意k柱上游|x_k|处(x_k<0)，把k柱以及k柱下游各柱的绕射势求和得到：

其中，

A_n为第1柱上的系数，i是虚数单位，ω是波浪圆频率，t是时间，n是整数，Z_n＝J_n′(Ka)/H_n′(Ka)，K为波数，J_n是第一类贝塞尔函数，H_n是第一类汉克尔函数，a是圆柱半径；

将汉克尔函数取渐近表达式进行简化后有：

其中，

上式第二项和第三项等价于两个与k柱左传波波程差为2(R/2)和2(N+1/2-k)R的等价圆柱左传波的叠加效果，这两个等价圆柱位于k柱下游R/2和(N+1/2-k)R处；

利用离k柱更远的等价圆柱的位置对初步波动间距表达式进行修正后即可得到波浪力曲线区域III中波动间距的最终表达式：

同样的，修正后的波浪力曲线区域III中任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式为：

在本发明的一个实施方式中，对于波浪入射角不等于零的情况，任意峰点横坐标的最终表达式一、任意谷点横坐标的最终表达式一和最终波动间距表达式一的获取过程如下：

其中，k柱为陈列中的任意一个圆柱，β为波浪入射角，令β≠0时，在区域III中任意第s和第s+1个峰点、谷点对应的柱间距是

将这些表示柱间距的符号再增加上角标u后，表示由上游各柱右转波作用所导致的相应的量；

首先得到任意峰点横坐标的初步达式一、任意谷点横坐标的初步表达式一和初步波动间距表达式一：

然后对上述三式进行修正即得到任意峰点横坐标的最终表达式一、任意谷点横坐标的最终表达式一和最终波动间距表达式一；

在本发明的一个实施方式中，对于波浪入射角不等于零的情况，任意峰点横坐标的最终表达式二、任意谷点横坐标的最终表达式二和最终波动间距表达式二的获取过程如下：

其中，

等表示柱间距的符号再增加上角标l表示由下游各柱左传波作用所导致的相应的量；

首先得到任意峰点横坐标的初步表达式二、任意谷点横坐标的初步表达式二和初步波动间距表达式二：

然后对上述三式进行修正，即得到任意峰点横坐标的最终表达式二、任意谷点横坐标的最终表达式二和最终波动间距表达式二：

在本发明的一个实施方式中，所述步骤600中的最终表达式在波浪入射角等于零时，任意圆柱波浪力曲线在区域III中波动间距为：

在波浪入射角不等于零时，任意圆柱波浪力曲线在区域III中波动间距为：

在本发明的一个实施方式中，所述区域III的自然数取值范围在2～10之间，所述区域II的自然数取值范围在5～10之间，所述区域I的自然数取值范围在40～50之间；当所述区域III的计算步长取五分之一的区域III最小波动间距，所述区域II的计算步长取十分之一的区域III最小波动间距，所述区域I的计算步长取五十分之一的区域III最小波动间距时，波浪力曲线的计算精度达到相对误差在1％以内。

在本发明的一个实施方式中，获取各所述峰点的横坐标的过程如下：

对于波浪入射角等于零的情况，波浪力曲线区域III中任意峰点横坐标的以柱间距-波长比表示的最终表达式为:

对于波浪入射角不等于零的情况，结合任意峰点横坐标的最终表达式一和任意峰点横坐标的最终表达式二可以得到以柱间距-波长比表示的任意峰点横坐标的最终表达式：

获取各所述谷点的横坐标的过程如下：

对于波浪入射角等于零的情况，波浪力曲线区域III中任意谷点横坐标的以柱间距-波长比表示的最终表达式为:

对于波浪入射角不等于零的情况，结合任意谷点横坐标的最终表达式一和任意谷点横坐标的最终表达式二可得到以柱间距-波长比表示的任意谷点横坐标的最终表达式：

在本发明的一个实施方式中，每一个所述峰点坐标和所述谷点坐标中的纵坐标获取方式如下：

根据坐底圆柱阵列水波绕射问题中任意k柱附近速度势的空间因子φ(r_k,θ_k)公式：

其中，

为绕射系数，k为圆柱阵列中任意一根圆柱的编号，编号k增大方向与x轴正方向保持一致，(r_k,θ_k)为垂直轴z轴通过k柱轴线的局部圆柱坐标系的极坐标，Z_n＝J_n′(Ka)/H_n′(Ka)，K为波数，a为圆柱半径，J_n为第一类贝塞尔函数，H_n是第一类汉克尔函数，n为整数。

求解所述速度势表达式中绕射系数的线性方程组如下：

其中，β为平面入射波传播方向与圆柱阵列整体坐标系中x轴正方向的夹角(波浪入射角)，且整体坐标系的建立使波浪入射角β≤π/2，R_jk为第k个柱子到第j个柱子的距离，i为虚数单位，m为整数,α_jk为第k个柱子到第j个柱子的方向角，I_k为入射波在第k个柱子的相位因子；

将所述峰点和所述谷点对应的波数K代入该方程可求得所述峰点和所述谷点对应波数下绕射系数

的值，将绕射系数

代入下式

就可得到波浪力曲线在区域III中任意峰点和谷点对应波数下，任意圆柱k所受沿圆柱阵列水平截面中各圆心连线方向的波浪力F^k；其中，ρ为水的密度，g为重力加速度，A为入射波的波幅，h为水深；

用相同几何尺寸的圆柱在相同环境条件下所受波浪力对前式所示的波浪力进行无量纲化，可得到波浪力曲线在区域III中任意峰点和谷点对应波数下圆柱阵列中任意第k柱的无量纲波浪力为：

对该式取模即得到无量纲波浪力幅值，这也就是波浪力曲线在区域III中任意峰点和谷点处的纵坐标。

在本发明的一个实施方式中，所述步骤200至步骤600中描述模型的获取过程还可采用如下获取方法：

步骤201，建立由排成一条直线的有限个相同坐底圆柱所组成的圆柱阵列的整体坐标系，列出坐底圆柱阵列水波绕射问题中任意圆柱附近的速度势表达式，速度势表达式中由线性方程组确定的未知系数称为绕射系数，由于绕射系数幅值随无量纲波数变化形成的绕射系数曲线与波浪力幅值随无量纲波数变化形成的波浪力曲线的波动特性和波动间距一致，因此对波浪力曲线的研究转化为对绕射系数曲线的研究；针对确定绕射系数的线性方程组采用克莱姆法则求解，用上述线性方程组右端列向量替换该线性方程组左端系数方阵相应的列然后得到目标方阵行列式；

步骤301，忽略高阶小量得到目标方阵行列式模的平方表达式；忽略平方表达式中随无量纲波数缓慢变化的慢变函数，得到与波动间距分析密切相关的目标表达式；对目标表达式求一阶导数，令该一阶导数为零，得到超越方程；

步骤401，利用超越方程分别推导波浪入射角等于零和不等于零时任意圆柱波浪力曲线在区域III中的波动间距表达式，综合各表达式后即得到区域III的描述模型。

在本发明的一个实施方式中，所述步骤201中的速度势表达式如下：

其中，

为绕射系数，k为圆柱阵列中任意一根圆柱的编号，编号k增大方向与x轴正方向保持一致，(r_k,θ_k)为垂直轴z轴通过k柱轴线的局部圆柱坐标系的极坐标，Z_n＝J_n′(Ka)/H_n′(Ka)，K为波数，a为圆柱半径，J_n为第一类贝塞尔函数，H_n是第一类汉克尔函数，n为整数。

在本发明的一个实施方式中，求解所述速度势表达式中绕射系数的线性方程组如下：

其中，β为波浪入射角，是由平面入射波传播方向与圆柱阵列整体坐标系中x轴正方向形成的夹角，且整体坐标系的建立使波浪入射角β≤π/2，K为波数,R_jk为第k个柱子轴线到第j个柱子轴线的距离，i为虚数单位，m为整数,N为圆柱阵列中的圆柱总数，α_jk为第k个柱子到第j个柱子的方向角，I_k为入射波在第k个柱子的相位因子；

利用克莱姆法则求解上述线性方程组后得到的绕射系数如下：

其中，D为所述线性方程组中系数方阵的行列式，

为用所述线性方程组的右端列向量替换系数方阵对应列得到的目标方阵行列式。

在本发明的一个实施方式中，对绕射系数

的分析可进一步转化为对目标方阵行列式模的平方

的分析，将

展开并忽略高阶小量，进而求其模的平方后得到的公式即为目标方阵行列式模的平方表达式：

其中，

κ＝Kd/π为无量纲波数，2d为相邻圆柱轴线之间的距离，M为截断项数。

在本发明的一个实施方式中，所述目标表达式为：

对目标表达式求一阶导数，令该一阶导数为零，得到的所述超越方程为：

[2(N-k)+1]tan(2πκ)＝tan{[2(N-k)+1]2πκ}。

在本发明的一个实施方式中，利用超越方程推导任意圆柱k在波浪入射角等于零时其在区域III中的波浪力曲线波动间距的表达式过程为：

步骤410，对超越方程左右两端表达式分别做泰勒展开以得到各自的近似表达式；

步骤411，再把各近似表达式分别代入超越方程，解析求得任意圆柱k在区域III波浪力曲线波动间距中的上下限表达式；

步骤412，忽略小量后，上述上下限表达式相同，因而可得到波浪入射角等于零时，圆柱阵列中任意圆柱k波浪力曲线在区域III中的波动间距描述模型的表达式为

在本发明的一个实施方式中，利用超越方程推导任意圆柱k在波浪入射角不等于零时波浪力曲线在区域III中的波动间距表达式时，需要将超越方程[2(N-k)+1]tan(2πκ)＝tan{[2(N-k)+1]2πκ}中的2πκ用(1+cosβ)πκ、(1-cosβ)πκ代替以得到波浪入射角不等于零时的两个修正超越方程，然后对这两个修正超越方程的左、右两端表达式分别做泰勒展开得到各自近似表达式；再把近似表达式分别代入两个修正超越方程，则可解析求得区域III波浪力曲线波动间距上下限的表达式，进而可以得到波浪入射角不等于零情况下的波浪力曲线波动间距表达式。

在本发明的一个实施方式中，在波浪入射角不等于零时，波浪力曲线按波动特性以及波动间距不同分为以下五种情形：

(1)[2(N-k)+1](1+cosβ)＞＞[2(k-1)+1](1-cosβ)；

(2)[2(k-1)+1](1-cosβ)＞＞[2(N-k)+1](1+cosβ)；

(3)[2(N-k)+1](1+cosβ)>[2(k-1)+1](1-cosβ)＞＞4；

(4)[2(k-1)+1](1-cosβ)>[2(N-k)+1](1+cosβ)＞＞4；

(5)[2(N-k)+1](1+cosβ)和[2(k-1)+1](1-cosβ)量级比较接近，均远大于4；

其中的“＞＞4”和“远大于4”表示区域III的波动间距远小于区域I各个高耸尖峰之间的距离，

第(1)种情形的波浪力曲线区域III波动间距表达式为：

第(2)种情形的波浪力曲线区域III波动间距表达式为：

第(3)种情形的波浪力曲线区域III波动间距表达式采用第(1)种情形和第(2)种情形间距表达式中的波动间距更小者；

第(4)种情形的波浪力曲线区域III波动间距表达式采用第(1)种情形和第(2)种情形间距表达式中的波动间距更小者；

第(5)种情形的波浪力曲线区域III波动间距表达式采用第(1)种情形和第(2)种情形间距表达式中的波动间距更小者。

在本发明的一个实施方式中，综合在波浪入射角不等于零时五种情形下区域III的波动间距表达式，得到波浪入射角不等于零时圆柱阵列中任意圆柱k波浪力曲线在区域III中最小波动间距

描述模型的表达式为：

本发明提供了这样的认识和理解：波浪力幅值随无量纲波数变化形成的波浪力曲线中非捕获区(区域III)的波动间距不随无量纲波数的变化而改变，该波动间距只与圆柱阵列中圆柱总数、标识圆柱位置的柱子编号和波浪入射角有关，可以用本发明给出的公式准确预测。

本发明可以加深对波浪力幅值随无量纲波数变化形成的波浪力曲线中非捕获区(区域III)波动特性的理解和认识，本发明提供了非捕获区中波动间距的预测公式。基于本发明，在进行有关工程结构设计和评估时，可以在保证精度的前提下减少工作量，缩短设计和评估周期，为提高工程结构的设计和评估水平提供技术支撑。

采用本发明给出的波浪力曲线计算步长确定方法，只要已知圆柱阵列中圆柱总数、标识圆柱位置的圆柱编号和波浪入射角，就可以在任意圆柱波浪力曲线计算开始之前，准确预测出计算步长取值多大即可得到准确的波浪力曲线。从而可以在不盲目增加不必要计算时间的前提下，得到准确的波浪力曲线，也就可以在保证精度的前提下减少工作量，缩短设计和评估周期，为提高工程结构的设计和评估水平提供技术支撑。

本发明给出了圆柱阵列波浪力曲线区域III中上下包络线的确定方法，采用本发明给出的峰点谷点横纵坐标的计算方法，可以直接预先确定出波浪力曲线区域III中峰点和谷点的位置，得到波浪力曲线区域III中的上下包络线。因而在方案评估、选型或初步设计时，可以在不漏掉波浪力曲线峰谷点的同时，避免计算对实际工程的这几个阶段来说不太重要的大量其他数据点，从而可以大幅减少计算量，节约成本、提高效率，为最终提高评估和设计水平提供技术支撑。

附图说明

图1为本发明一个实施方式中相同直径圆柱排成一条直线所组成的圆柱阵列示意图；

图2为圆柱总数N＝17，柱子编号k＝9，波浪入射角β＝0，直径-柱间距比a/d＝1/4的单排坐底柱群波浪力幅值随无量纲波数变化形成的波浪力曲线和捕获有关区域(RegionI和Region II)和非捕获区(Region III)等三个区域的示意图；

图3为圆柱总数N＝301，柱子编号k＝151，波浪入射角β＝0，直径-柱间距比a/d＝1/2的单排坐底圆柱阵列在near-trapping尖峰对应的无量纲波数下波浪力幅值随柱子编号变化图像；

图4为圆柱总数N＝301，柱子编号k＝151，波浪入射角β＝0，直径-柱间距比a/d＝1/2的单排坐底圆柱阵列在near-trapping尖峰左边第一个谷点对应的无量纲波数下的波浪力幅值随柱子编号变化图像；

图5为圆柱总数N＝301，柱子编号k＝151，波浪入射角β＝0，直径-柱间距比a/d＝1/2的单排坐底圆柱阵列在near-trapping尖峰左边第一个峰点对应的无量纲波数下的波浪力幅值随柱子编号变化图像；

图6为波浪相长干涉/相消干涉的示意图；其中，波动曲线为N＝21，柱子编号k＝1，波浪入射角β＝0，直径-柱间距比a/d＝1/4的单排坐底圆柱阵列波浪力幅值随无量纲波数的变化曲线；

图7为圆柱总数N＝101，波浪入射角β＝0，直径-柱间距比a/d＝1/4，柱子编号k不同的波浪力幅值波动间距测量值随着无量纲波数变化的图；

图8为圆柱数目为N＝11,21,51,101，波浪入射角β＝0，直径柱间距比a/d＝1/4时，第k柱受到的波浪力幅值在区域III的波动间距测量值

和理论计算值

随柱子编号变化图像；

图9为本发明另一实施方式的描述模型流程示意图；

图10为γ(κ)和χ(κ)以及它们的泰勒展开式的多个交点的示意图；

图11为圆柱总数N＝301，直径-柱间距比a/d＝1/4，波浪入射角不等于0时，波动间距的五种情形的波浪力幅值随无量纲波数变化形成的波浪力曲线(左列)以及快速傅里叶变换得到的数值结果(右列)。

具体实施方式

在详细阐述本发明具体内容之前，先介绍综合现有研究成果和我们的分析结果得到的确定波浪力曲线捕获有关区域(区域I和区域II)位置和范围的方法。(区域I和区域II即为前文所述捕获有关区域，这两个区域具有明显比区域III更高更深的峰谷，此外，这两个区域的另一个特征是：波动间距随无量纲波数改变而发生变化)

1)区域I(捕获有关区域)

目前已有不少文献研究了无限长圆柱阵列或水槽中心线上布置单个圆柱的trappedmode频率问题，这些结果可以估算有限数目坐底圆柱阵列的near-trapping波数，即可以得到区域I的位置。具体做法是，根据圆柱直径-柱间距之比a/d(2a为圆柱直径，2d为相邻圆柱轴线之间的距离)，查找文献中已知的trapped mode对应的波数，在该波数附近搜索计算找到尖峰便可以得到有限长圆柱阵列的捕获有关区域中的区域I。对于某些a/d文献中未给出对应trapped mode波数的情况，可以采用1/[20(N-k)+10]作为一个初始计算步长(N为圆柱阵列中圆柱总数，k为标识柱子位置的柱子编号)，在Kd/π为0.5的整数倍附近搜索计算(K为波数)找到尖峰便可以得到有限长圆柱阵列的捕获有关区域中的区域I。对于有限长单排圆柱阵列，随着直径-柱间距比a/d减小，区域I的尖峰点也向右移动。通过与接近的a/d所对应的结果对比，可以进一步缩小区域I的计算范围。对于所得到尖峰对应的波数值，绘制出波浪力幅值与柱子编号的关系图，如能呈现一个完整的半波形式，如图3所示，最大的波浪力作用在中间柱上，则说明该尖峰位置足够准确。如果不是如此，还可以继续加密步长找到更准确的尖峰点。

2)区域II(另一个捕获有关区域)

在波浪力曲线区域I附近的次级峰谷构成了区域II，区域II中曲线的波动间距随无量纲波数改变而发生变化。目前已有文献对有限长圆柱阵列波浪力曲线中区域I尖峰左边的一些次级峰、谷进行了研究。研究表明，这些次级峰、谷与无限长圆柱阵列Rayleigh-Bloch波问题以及水槽中横向排布多柱体的Trapped modes有关。特别地，对于由N个单排坐底圆柱组成的圆柱阵列的中间柱，其波浪力曲线尖峰点左边的次峰、次谷位置的横坐标(无量纲波数)与圆柱数目为N/2，N/3，N/4…的圆柱阵列之中间柱波浪力曲线尖峰位置的横坐标(无量纲波数)严格对应，具体如下：

圆柱数目为N/2的单排坐底圆柱阵列之中间柱波浪力曲线尖峰位置横坐标对应着圆柱数目为N的圆柱阵列之中间柱波浪力曲线尖峰左边第一个谷点位置的横坐标，此谷点位置横坐标对应的无量纲波数下圆柱数目为N的圆柱阵列的波浪力幅值与圆柱编号关系图呈现出两个半波的形式，如图4所示，这两个半波的最高峰对应的波浪力幅值与同样波数下圆柱数目为N/2的圆柱阵列之中间柱的波浪力幅值大小相等。

圆柱数目为N/3的单排坐底圆柱阵列之中间柱波浪力曲线尖峰位置横坐标对应着圆柱数目为N的圆柱阵列之中间柱波浪力曲线尖峰左边第一个峰点位置的横坐标，此峰点位置横坐标对应的无量纲波数下圆柱数目为N的圆柱阵列的波浪力幅值与圆柱编号关系图呈现出三个半波的形式，如图5所示，这三个半波的最高峰对应的波浪力幅值与同样波数下圆柱数目为N/3的圆柱阵列之中间柱的波浪力幅值大小相等。

圆柱数目为N/4，N/5…等的情况与上述类似，可依此类推。通常来说，当N/n_i～10时(n_i为自然数)，near-trapping的影响已经相当微弱，可将此时圆柱数目为(N/n_i＝)10的圆柱阵列波浪力曲线尖峰位置对应的无量纲波数作为N个圆柱组成圆柱阵列波浪力曲线区域II的左边界限。

我们通过计算分析发现，对于不同的直径-柱间距比a/d，区域II受到near-trapping的影响范围不同。a/d越大，near-trapping的影响范围亦越大。例如，对于a/d＝0.25的情况，圆柱数目N/n_i～20阵列的波浪力尖峰位置对应的无量纲波数可作为区域II的左界限，而对于a/d＝0.5的情形，这个左界限会持续到圆柱数目N/n_i～5阵列的波浪力尖峰位置对应的无量纲波数。对于单排坐底圆柱群阵列中任意第k柱的情况，可以参照上述中间柱的范围来确定。

本发明中的圆柱阵列是指，较大数量穿透水面的相同直径圆柱排成一条直线所组成的圆柱阵列(即，该圆柱阵列的水平截面中各个圆心在一条直线上)。本发明中的波浪力是指，任意圆柱所受沿圆柱阵列水平截面中各圆心连线方向的波浪力。本发明中的波动间距是指，波浪力幅值随无量纲波数变化形成的波浪力曲线上，相邻两个极大值点(或极小值点)横坐标之间的距离。在本发明中，也用“峰”或“谷”来描述极大值点或极小值点。

本发明提供圆柱阵列波浪力幅值波动特性、峰谷和包络线的获取方法，包括如下步骤：

步骤100，将波浪力幅值随无量纲波数变化形成的波浪力曲线中的多个高耸尖峰所在区域作为区域I，将高耸尖峰附近低于所述高耸尖峰且曲线波动间距随无量纲波数改变而变化的次级峰和谷所在区域作为区域II，将排除区域I和区域II后的波浪力曲线作为区域III；

区域III(Region III)称为非捕获区，区域I(Region I)和区域II(Region II)称为捕获有关区域，在非捕获区(区域III)，波浪力曲线具有非常有规律的波动现象。

步骤200，建立由排成一条直线的有限个相同坐底圆柱所组成的圆柱阵列坐标系并确定相关参数，将无量纲波数转化为相邻两圆柱间距离和波长的比值，则各圆柱的绕射波发生相长干涉时波浪力曲线出现峰点，发生相消干涉时波浪力曲线出现谷点；对于波浪入射角等于零的情形，分析入射波传到任意某个圆柱发生绕射及入射波传到该圆柱下游最后一个圆柱发生绕射这两个路径的波程差，求出波浪力曲线区域III中任意峰点横坐标的初步表达式和相邻峰点横坐标之差的初步表达式，以及任意谷点横坐标的初步表达式和相邻谷点横坐标之差的初步表达式，进而可知相邻两个峰点或相邻两个谷点的横坐标之差相等，最后得到任意圆柱波浪力曲线在区域III中的初步波动间距表达式；

波程差与相邻两圆柱间距离成正比，如图6所示，令该波程差分别等于s倍和s+1倍的波长(s为任意自然数)，可求得波浪力曲线区域III中任意峰点横坐标的初步表达式和相邻峰点横坐标之差的初步表达式，令该波程差分别等于2s-1倍和2s+1倍半波长，可求得波浪力曲线区域III中任意谷点横坐标的初步表达式以及相邻谷点横坐标之差的初步表达式；

这里圆柱阵列中圆柱的数量通常大于9个，波浪力幅值曲线的相邻极大或相邻极小值点的横坐标间距在区域III中是定值，不随无量纲波浪频率改变，只与阵列中圆柱总数N、标识圆柱位置的柱子编号k、波浪入射角β有关，而且可以非常准确地用简单公式来预测。

在圆柱阵列坐标系中，平面入射波传播方向与圆柱阵列整体坐标系中x轴正方向的夹角称为波浪入射角，且整体坐标系的建立使波浪入射角小于等于90度；k为圆柱阵列中的任意一根圆柱的编号，编号k增大方向与圆柱阵列整体坐标系中x轴正方向保持一致。

其中的初步表达式获取过程如下：

无量纲波数(即波浪力曲线横坐标)Kd/π可以改写为：Kd/π＝R/λ；于是，波浪力幅值随无量纲波数变化也可以认为是波浪力幅值随着柱间距—波长比而变化。当j柱的绕射波传到k柱附近与k柱绕射波发生相长(相消)干涉时，k柱上的波浪力幅值会得到一个峰(谷)值。R_p(1)和R_v(1)分别表示波浪力曲线中的第一个峰点和第一个谷点对应的柱间距，对于波浪入射角β＝0的情况，波浪力曲线中的第一个峰点R_p(1)和第一个谷点R_v(1)对应的柱间距用下式表示：

2(N-k)R_p(1)＝λ

其中，K是波数，R＝2d是相邻圆柱的柱间距，λ是波长，脚标中，p代表峰点，v代表谷点，(1)代表第一个峰或者谷。

令在区域III中任意第s和第s+1个峰点对应的柱间距是R_p(s)和R_p(s+1)，R_p(s+1)＝R_p(s)+δR_p，峰点出现的条件是发生相长干涉，则相邻的第s和第s+1个峰点对应的波程差应该分别是s倍和s+1倍的波长，于是N柱与k柱左传绕射波(以下简称“左传波”)的波程差应满足：

2(N-k)R_p(s)＝sλ

2(N-k)(R_p(s)+δR_p)＝(s+1)λ

上述两式相减得到相邻峰点横坐标之差的初步表达式为：

谷点出现的条件是发生相消干涉，对于第s和第s+1个谷点对应的柱间距R_v(s)和R_v(s+1)＝R_v(s)+δR_v，相应的波程差分别等于2s-1倍和2s+1倍半波长，因此有：

得到相邻谷点横坐标之差的初步表达式为：

由上述推导可知，相邻峰点和相邻谷点的横坐标之差相等，则初步波动间距表达式为：

N为圆柱阵列中的圆柱总数；此时波浪力曲线区域III中任意峰点横坐标的初步表达式和任意谷点横坐标的初步表达式为：

步骤300，在任意某个圆柱上游处，将该圆柱产生的向上游传播的左传绕射波和该圆柱下游各个圆柱产生的左传绕射波求和，利用汉克尔函数取渐近表达式进行简化，得到对波动间距问题来说与当前圆柱阵列等价、且位置与现有圆柱位置均不相同的两个位于该圆柱下游的等价圆柱的具体位置，将离该圆柱更远的那个等价圆柱的位置代入初步波动间距表达式、任意峰点横坐标的初步表达式和任意谷点横坐标的初步表达式进行修正，即可得到波浪入射角等于零时波浪力曲线区域III中波动间距的最终表达式、任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式；

获取波浪入射角等于零时波浪力曲线区域III中波动间距的最终表达式的过程如下：

在任意k柱上游|x_k|处(x_k<0)，把k柱以及k柱下游各柱的绕射势求和(就相当于前面所述左传波求和)得到：

其中，对于一长排柱子远离两端的中间区域柱子，有

A_n为第1柱上的系数，i是虚数单位，ω是波浪圆频率，t是时间，n是整数，Z_n＝J_n′(Ka)/H_n′(Ka)，K为波数，J_n是第一类贝塞尔函数，H_n是第一类汉克尔函数，a是圆柱半径；

将汉克尔函数取渐近表达式对前式进行简化后有：

该式中第二、三项的常数±π/2对于下面讨论的波动间距问题没有贡献；其中，

上式实际上就是三个左传平面波的叠加，第一项是k柱左传波，第二项和第三项是k+1柱到N柱这N-k个柱子的左传波之和，第二项和第三项等价于两个与k柱左传波波程差为2(R/2)和2(N+1/2-k)R的等价圆柱左传波的叠加效果，这两个等价圆柱位于k柱下游R/2和(N+1/2-k)R处；

利用离k柱更远的等价圆柱的位置对初步波动间距表达式进行修正后即可得到波浪力曲线区域III中波动间距的最终表达式：

即：

同样的，在入射角等于零时，修正后的波浪力曲线区域III中任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式为：

步骤400，对于波浪入射角不等于零的情形，利用前述同样方法，首先分析入射波传到任意某个圆柱发生绕射及入射波传到该圆柱上游的阵列端部第一个圆柱发生绕射这两个路径的波程差，然后分析入射波传到任意某个圆柱发生绕射及入射波传到该圆柱下游的最后一个圆柱发生绕射这两个路径的波程差，即可分别得到任意圆柱波浪力曲线在区域III中由于该圆柱上游圆柱绕射波的作用导致的任意峰点横坐标的初步表达式一、任意谷点横坐标的初步表达式一和初步波动间距表达式一，以及由于该圆柱下游圆柱绕射波的作用导致的任意峰点横坐标的初步表达式二、任意谷点横坐标的初步表达式二和初步波动间距表达式二；

入射角不为零(β≠0)的情况要比入射角等于零的情况复杂，此时k柱上游各柱右传绕射波与k柱右传绕射波(以下简称“右传波”)有了波程差，这一点与β＝0的情况不同。令β≠0时，在Region III中任意第s和第s+1个峰、谷点对应的柱间距是

由于需要分别与考虑上游各柱右传波和下游各柱左传波有关的情况，在下文中，上述表示柱间距的符号当增加上角标“u”时表示由上游各柱右传波作用所导致的相应的量，当增加上角标“l”时表示由下游各柱左传波作用所导致的相应的量。

对于波浪入射角不等于零的情况，首先分析位于k柱上游的阵列端部第一个圆柱向下游传播的右传绕射波与k柱右传绕射波的波程差，使其满足发生相长或相消干涉的条件，可得到初步波动间距表达式一，过程如下：

其中，k柱为圆柱阵列中的任意一个圆柱，β为波浪入射角，得到初步波动间距表达式一：

β≠0时，Region III中任意第s个峰点横坐标的初步表达式一和任意第s个谷点横坐标的初步表达式一为：

再分析位于k柱下游的最后一个圆柱向上游传播的左传绕射波与k柱左传绕射波的波程差，使其满足发生相长或相消干涉的条件，可得到初步波动间距表达式二，过程如下：

其中，l表示与下游各柱左传波有关的量；得到初步波动间距表达式二如下：

Region III中任意第s个峰点横坐标的初步表达式二和任意第s个谷点横坐标的初步表达式二为：

步骤500，首先在任意某个圆柱的下游处，对该圆柱产生的向下游传播的右传绕射波和该圆柱上游各个圆柱产生的右传绕射波求和，利用前述同样的方法修正，得到波浪入射角不等于零时，任意峰点横坐标的最终表达式一、任意谷点横坐标的最终表达式一和最终波动间距表达式一；然后，在任意某个圆柱的上游处，对该圆柱产生的向上游传播的左传绕射波和该圆柱下游各个圆柱产生的左传绕射波求和，同样方法修正后得到任意峰点横坐标的最终表达式二、任意谷点横坐标的最终表达式二和最终波动间距表达式二；

与β＝0过程类似，在k柱下游|x_k|处(x_k>0)，k柱以及k柱上游各柱的绕射势之和为

其中，

该式中，第二项和第三项表示第1柱到k-1柱这k-1个柱子的右传波之和。这两项的效果可以看作是位于第k柱上游R/2和(k-1/2)R处的两个柱子的效果，于是，初步波动间距表达式一

经过修正后得到入射角不等于零时最终波动间距表达式一

及此状态下的任意峰点横坐标的最终表达式一和任意谷点横坐标的最终表达式一为；

在k柱上游|x_k|处(x_k<0)，k柱以及k柱下游各柱的绕射势之和为：

其中，

该式中，第二项和第三项表示第k+1柱到N柱这N-k个柱子的左传波之和，这两项的效果可以看作是位于k柱下游R/2和(N+1/2-k)R处的两个圆柱的左传波的效果，于是波浪入射角不等于零时的初步间距表达式二

经过修正后得到最终波动间距表达式二：

及此状态下的任意峰点横坐标的最终表达式二和任意谷点横坐标的最终表达式二为；

步骤600，采用最终波动间距表达式一和最终波动间距表达式二给出的波动间距的更小者，就得到波浪入射角不等于零时任意圆柱波浪力曲线在区域III中最小波动间距的最终表达式；综合波浪入射角等于零时和不等于零时的波动间距的最终表达式，即可得到任意圆柱波浪力曲线在区域III中波动间距描述模型的最终表达式；

最终表达式在波浪入射角等于零时，任意圆柱波浪力曲线在区域III中波动间距为：

在波浪入射角不等于零时，任意圆柱波浪力曲线在区域III中波动间距为：

图7是柱子总数N＝101,波浪入射角β＝0，直径柱间距比a/d＝1/4时，波浪力幅值的波动间距测量值

随着无量纲波数变化的曲线。可以看到，波动间距在很大波数范围下是恒定的，这个区域就是区域III，之后波动间距迅速下降的区域就是区域I和II，图中渐近线及数值是利用本发明最终表达式计算得到的理论预测值，结果非常吻合。

图8分别是N＝11,N＝21,N＝51,N＝101情况下，波浪入射角β＝0，直径柱间距比a/d＝1/4时，单排坐底圆柱群阵列第k柱受到的波浪力在区域III的波动间距测量值

和描述模型表达式计算值

的对比。通过对比发现，本发明的描述模型表达式预测值与实际计算值吻合非常好。

步骤700，根据入射角等于零和不等于零时波动间距的最终表达式计算出区域III中的最小波动间距，以该最小波动间距作为波浪力曲线在区域III中计算步长的上限，根据不同精度要求将最小波动间距除以相应的自然数作为区域III中计算步长的下限，从而得到区域III的计算步长；

根据所需精度要求将区域III的最小波动间距除以相应精度要求对应的自然数，得到区域II中的计算步长；根据所需精度要求将区域III的最小波动间距除以相应精度要求对应的自然数，得到区域I中的计算步长；

在入射角等于零时波动间距的最终表达式

和

入射角不等于零时波动间距的最终表达式

中代入N，k和β后即可计算出任意圆柱波浪力曲线在区域III中的最小波动间距

或

以下为简洁起见，统一用符号

表示波浪入射角等于零和不等于零两种情况下的最小波动间距，需注意的是波浪入射角β＝0时采用

的表达式计算最小波动间距；此时的最小波动间距

作为区域III计算步长的上限，而计算步长的下限，则根据计算时的精度要求选取一个范围在2～10之间的自然数，作为最小波动间距

的除数，即可确定该下限，这里自然数的取值越大，则精度越高，相应花费的计算时间也越长。通过确定后的上限和下限，即可得到区域III的计算步长。

由于区域II范围比区域III范围小很多且波动间距变小，因此，可根据精度要求取一个范围在5～10之间的自然数作为区域III的最小波动间距

的除数，即可得到区域II的计算步长，这里同样是自然数的取值越大，则精度越高，相应花费的计算时间也越长。

由于区域I范围比区域II小很多且具有分离的高耸尖峰，因此，可根据精度要求取一个范围在40～50之间的自然数作为区域III的最小波动间距

的除数，即可得到区域I的计算步长，这里同样是自然数的取值越大，则精度越高，相应花费的计算时间也越长。

具体地，当区域III的计算步长取五分之一的区域III最小波动间距

区域II的计算步长取十分之一的区域III最小波动间距

区域I的计算步长取五十分之一的区域III最小波动间距

时，波浪力曲线的计算精度达到相对误差在1％以内。

步骤800，采用最终表达式一和最终表达式二中横坐标的更小者作为波浪入射角不等于零时任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式，综合波浪入射角等于零时和不等于零时的最终表达式，即可得到任意圆柱波浪力曲线在区域III中任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式；根据横坐标的最终表达式得到任意峰点和谷点处的横坐标，从而可得对应的波数，求解线性方程组以获得速度势表达式中的未知绕射系数，进而可得到任意圆柱所受波浪力，将该波浪力无量纲化并取模，就得到波浪力曲线在区域III中任意峰点和谷点处的纵坐标；

根据上述获取的任意峰点和谷点处的横坐标以及任意峰点和谷点处的纵坐标，确定出每一个峰点和谷点的位置，将这些峰点连线，就得到波浪力曲线在区域III中的上包络线；将这些谷点连线，就得到波浪力曲线在区域III中的下包络线。

任意谷点横坐标的最终表达式一和最终表达式二分别体现了上游柱子和下游柱子对谷点横坐标的影响，

结合此二式可得到以柱间距-波长比表示的任意谷点横坐标的最终表达式：

任意峰点横坐标的最终表达式一和最终表达式二分别体现了上游柱子和下游柱子对峰点横坐标的影响，

结合此二式可得到以柱间距-波长比表示的任意峰点横坐标的最终表达式：

每一个峰点坐标和谷点坐标中的纵坐标获取方式如下：

计算坐底圆柱阵列水波绕射问题中任意k柱附近速度势的空间因子φ(r_k,θ_k)的公式为：

其中，

为绕射系数，k为圆柱阵列中任意一根圆柱的编号，编号k增大方向与x轴正方向保持一致，(r_k,θ_k)为垂直轴z轴通过k柱轴线的局部圆柱坐标系的极坐标，Z_n＝J_n′(Ka)/H_n′(Ka)，K为波数，a为圆柱半径，J_n为第一类贝塞尔函数，H_n是第一类汉克尔函数，n为整数。

其中的绕射系数

(未知系数)由以下绕射系数线性方程确定：

其中，β为平面入射波传播方向与圆柱阵列整体坐标系中x轴正方向的夹角(波浪入射角)，且整体坐标系的建立使波浪入射角β≤π/2，R_jk为第k个柱子到第j个柱子的距离，i为虚数单位，n,m为傅里叶模态截断项,α_jk为第k个柱子到第j个柱子的方向角，I_k为入射波在第k个柱子的相位因子。

由于横坐标R/λ等于无量纲波数Kd/π，可知任意谷点横坐标的最终表达式和任意峰点横坐标的最终表达式也就是以无量纲波数表示的峰点和谷点的横坐标，因而可以得到峰点和谷点对应的波数K。

将峰点和谷点对应的波数K代入前述绕射系数线性方程可求得峰点和谷点对应波数下，绕射系数

的值，将

代入下式

就可得到波浪力曲线在区域III中任意峰点和谷点对应波数下，任意圆柱k所受沿圆柱阵列水平截面中各圆心连线方向的波浪力F^k；其中，ρ为水的密度，g为重力加速度，A为入射波的波幅，h为水深。

用相同几何尺寸的圆柱在相同环境条件下所受波浪力对前述公式所示的波浪力F^k进行无量纲化，可得到波浪力曲线在区域III中任意峰点和谷点对应波数下，圆柱群阵列中任意第k柱的无量纲波浪力为：

对该式取模就得到无量纲波浪力幅值，这也就是波浪力曲线在区域III中任意峰点和谷点处的纵坐标。

至此，波浪力曲线在区域III中任意峰点和谷点处的横坐标和纵坐标均已得到，将这些谷点连线，就得到波浪力曲线在区域III中的下包络线；将这些峰点连线，就得到波浪力曲线在区域III中的上包络线。

如图9所示，在本发明的一个实施例中，提供另一种建立圆柱阵列的波浪力幅值在非捕获区波动间距的描述模型的方法，该方法与前述利用相长相消理论建立描述模型的方法区别在于步骤200至步骤600的处理步骤不同，而求取波浪力曲线计算步长的方式不变，以下为了方便理解，并未去掉与前述方法相同的步骤100，采用完整描述的方式，具体步骤如下：

步骤101，将波浪力幅值随无量纲波数变化形成的波浪力曲线中的多个高耸尖峰所在区域作为区域I，将高耸尖峰附近低于所述高耸尖峰且曲线波动间距随无量纲波数改变而变化的次级峰和谷所在区域作为区域II，将排除区域I和区域II后的波浪力曲线作为区域III进行分析；

区域III(Region III)称为非捕获区，区域I(Region I)和区域II(Region II)称为捕获有关区域，在非捕获区(区域III)，波浪力曲线具有非常有规律的波动现象。

步骤201，建立由排成一条直线的有限个相同坐底圆柱所组成的圆柱阵列的整体坐标系，列出坐底圆柱阵列水波绕射问题中任意圆柱附近的速度势表达式，速度势表达式中由线性方程组确定的未知系数称为绕射系数，由于绕射系数幅值随无量纲波数变化形成的绕射系数曲线与波浪力幅值随无量纲波数变化形成的波浪力曲线的波动特性和波动间距一致，因此对波浪力曲线的研究转化为对绕射系数曲线的研究；针对确定绕射系数的线性方程组采用克莱姆法则(Cramer's Rule)求解，用上述线性方程组右端列向量替换该线性方程组左端系数方阵相应的列然后得到目标方阵行列式；

这里圆柱阵列中圆柱的数量通常大于9个，波浪力幅值曲线的相邻极大或相邻极小值点的横坐标间距在区域III中是定值，不随无量纲波浪频率改变，只与阵列中圆柱总数N、标识圆柱位置的柱子编号k、波浪入射角β有关，而且可以非常准确地用简单公式来预测。

其中的速度势表达式如下：坐底圆柱阵列水波绕射问题中任意k柱附近速度势的空间因子φ(r_k,θ_k)为：

其中未知系数

称为绕射系数，由以下方程确定：

(1)、(2)式中，k为圆柱阵列中的任意一根圆柱的编号，编号k增大方向与圆柱阵列整体坐标系中x轴正方向保持一致，(r_k,θ_k)为垂直轴z轴通过k柱轴线的局部圆柱坐标系的极坐标，Z_n＝J_n′(Ka)/H_n′(Ka)，J_n为第一类贝塞尔函数，H_n是第一类汉克尔函数，n为整数；β为波浪入射角，是由平面入射波传播方向与圆柱阵列整体坐标系中x轴正方向形成的夹角，且整体坐标系的建立使波浪入射角β≤π/2，K为波数,a为圆柱半径，R_jk为第k个柱子轴线到第j个柱子轴线的距离，i为虚数单位，m为整数,N为圆柱阵列中的圆柱总数，α_jk为第k个柱子到第j个柱子的方向角，I_k为入射波在第k个柱子的相位因子。

开展有关计算和分析后，可以发现波浪力曲线区域III与绕射系数

幅值随无量纲波数变化形成的绕射系数曲线的区域III范围相同，波动间距也相同。因此对波浪力曲线波动间距的讨论可以转化为对绕射系数

随无量纲波数变化形成的绕射系数曲线的波动间距的研究。为了得到波动间距的描述模型，我们依据克莱姆法则(Cramer's Rule)求解绕射系数

线性方程组，绕射系数

可表示为

式中，D为式(2)线性方程组中系数方阵的行列式，

为用线性方程组的右端列向量替换系数方阵对应列得到的目标方阵行列式的值。

步骤301，忽略高阶小量得到目标方阵行列式模的平方表达式；忽略平方表达式中随无量纲波数缓慢变化的慢变函数，得到与波动间距分析密切相关的目标表达式；对目标表达式求一阶导数，令该一阶导数为零，得到超越方程；

经过研究发现，对绕射系数

的分析可进一步转化为对目标方阵行列式模的平方

的分析，将

展开并忽略高阶小量，进而求其模的平方后再忽略高阶小量，得到的公式即为目标方阵行列式模的平方表达式：

其中，

上述式中，κ＝Kd/π为无量纲波数，2d为相邻圆柱轴线之间的距离，M为截断项数。

对于波浪入射角β＝0的情况，根据式(4)-(6)可以得到

其中，

经过数值和理论分析表明，式(7)-(11)中，ξ(κ)，

和

是没有快速波动的慢变函数，对本发明讨论的波动间距不产生影响。α(κ)相对于随κ快速增长的4(j-k)πκ来说类似于一个常量，其存在只会引起函数的微小平移，对波动间距影响很小。于是，针对于式(7)的波动间距的研究可以用以下简化函数替代

此处，由于考察的是区域III的波动特性，则存在Kd/π≠μ/2(μ为整数)。即，在区域III中,sin2πκ≠0。根据三角公式，式(12)可以写成

其中，

上式中，随着j的增大，ε(j,κ)趋于零，可以忽略。于是式(13)进一步简化为

由于波动间距是相邻两个极大值点或者极小值点横坐标的水平距离，因此考察τ(κ)的一阶导数为零的条件，则得到超越方程

[2(N-k)+1]tan(2πκ)＝tan{[2(N-k)+1]2πκ} (16)

步骤401，利用超越方程分别推导波浪入射角等于零和不等于零时任意圆柱k波浪力曲线在区域III中的波动间距表达式，综合各表达式后即得到区域III的描述模型。

以下说明波浪入射角等于零(β＝0)和不等于零(β≠0)两种情况下的描述模型。

一、波浪入射角等于零时，单排坐底圆柱阵列波浪力曲线区域III中的波动规律如下：

为叙述方便，针对超越方程[2(N-k)+1]tan(2πκ)＝tan{[2(N-k)+1]2πκ}，定义γ(κ)＝[2(N-k)+1]tan(2πκ)和χ(κ)＝tan{[2(N-k)+1]2πκ}＝tan{[4(N-k)+2]πκ}。

式(16)对应着曲线γ(κ)和χ(κ)的交点，这两条曲线有多个交点，也就是说式(16)有多个解。由于[4(N-k)+2]π比2π大很多，因此在任一2πκ∈[pπ,(p+1)π]的区间，p＝0,1,2…，上述多个交点基本均在[4(N-k)+2]πκ＝qπ+π/2，q＝0,1,2…附近。图10(a)和(b)给出了交点的分布情况。显然，交点均在χ(κ)的渐近线附近。为了得到交点横坐标的解析表达式，首先把χ(κ)在qπ+π/2附近展开。于是在qπ+π/2，q＝0,1,2…附近有

把γ(κ)在pπ+π/2，p＝0,1,2…附近展开，有

把γ(κ)在在pπ附近展开有

如图10(b)和(c)所示，令γ(κ)和χ(κ)的交点为A，

和χ(κ)的交点为B，

和χ(κ)的交点为C。由于

交点A必然位于交点B以及交点C之间，即κ^B＞κ^A＞κ^C(或者κ^C＞κ^A＞κ^B)，κ^A,κ^B,κ^C分别为交点A,B,C的横坐标。得到κ^B,κ^C就得到了κ^A的上下界。实际上，绝大多数交点均在χ(κ)渐近线附近，由于渐近线附近导数很大，因而渐近线附近的点，即使纵坐标有明显区别，其横坐标的差别也非常小。

将式(17)和(18)代入式(16)，可以得到基于式(18)展开形式，式(16)的第q个解

将式(17)和(19)代入式(16)，可以得到基于式(19)展开形式，式(16)的第q个解

如图10(c)所示，对于κ^B＞κ^A＞κ^C的情况，有

根据式(20)，(21)和(22)，忽略小量，可以得到波浪入射角等于零时区域III波浪力曲线波动间距

的描述模型的表达式

二、波浪入射角不等零时单排坐底圆柱阵列波浪力曲线区域III中的波动规律如下：

对于波浪入射角β≠0的情况，波浪力幅值和未知系数

幅值的波动特性完全一样，此处仍然只讨论未知系数

幅值的波动特性。根据式(4)可得

其中，

υ＝2(j-k)πκcosβ+2|j-k|πκ (25)

上式中，将第k柱上游和下游的贡献分离，并利用式(13)–(15)的化简方式，有

其中，

与β＝0的情况不同的是，β≠0时，在不同参数组合下，区域III中波浪力曲线按波动特性可以分为以下五种情形：

情形1.[2(N-k)+1](1+cosβ)＞＞[2(k-1)+1](1-cosβ)

对于这种情形，

和

是慢变函数，而

和

是快变函数。此时，由第k柱下游柱子决定的

和

对最小波动间距起主要贡献，而

和

不影响最小波动间距。依照β＝0的情况类似的推导有，最小波动间距为：

下标“ds”表示第k柱下游柱子的贡献。

情形2.[2(k-1)+1](1-cosβ)＞＞[2(N-k)+1](1+cosβ)

与情形1相反，由第k柱上游柱子决定的

和

对最小波动间距起主要贡献，此时最小波动间距为

下标“us”表示第k柱上游柱子的贡献。

情形3.[2(N-k)+1」(1+cosβ)＞[2(k-1)+1」(1-cosβ)＞＞4

对于这种情形，尽管

和

比

和

振荡得慢，但是在相邻的两个区域I的尖峰之间能看出两种振荡频率的波动。较小的波动间距由式(31)确定，较大的波动间距由(32)确定。最小波动间距主要由来自下游柱子的贡献决定，也就是(31)和(32)式两者之间的更小者。

情形4.[2(k-1)+1](1-cosβ)＞[2(N-k)+1](1+cosβ)＞＞4

这种情形与情形3类似，最小波动间距主要由来自上游柱子的贡献决定，也就是(31)和(32)式两者之间的更小者。

情形5.[2(N-k)+1](1+cosβ)和[2(k-1)+1](1-cosβ)量级比较接近，均远大于4

这种情形下，区域III中波浪力曲线的波动性不是很规则，直观看起来似乎极值点间距失去了前面说的规律。但根据快速傅里叶变换的结果，实际上上述两种波动都存在，其内在机制没有改变，上述两种波动间距表达式同样采用式(31)和(32)，最小波动间采用(31)和(32)式两者之间的更小者。

上面所述“＞＞4”和“远大于4”表示区域III的波动间距远小于区域I各个高耸尖峰的间距。

综合以上五种情形，对于波浪斜射的情形(即波浪入射角不等于零的情况β≠0)，波浪力曲线最小波动间距

描述模型的表达式为：

图11给出了上述五种情形的例子，其中左侧一列是波浪力幅值随无量纲波数变化形成的波浪力曲线，图中给出了测量的结果和用描述模型表达式(33)计算得到的结果，吻合良好。作为进一步的验证，右侧一列给出了快速傅里叶变换给出的数值结果，快速傅里叶变换给出的频率和波动间距数值结果与描述模型表达式(33)的计算结果吻合良好。

综合波浪入射角等于零和不等于零时的两个表达式(23)和(33)即可得到波浪力曲线在区域III中的描述模型。

至此，本领域技术人员应认识到，虽然本文已详尽示出和描述了本发明的多个示例性实施例，但是，在不脱离本发明精神和范围的情况下，仍可根据本发明公开的内容直接确定或推导出符合本发明原理的许多其他变型或修改。因此，本发明的范围应被理解和认定为覆盖了所有这些其他变型或修改。

Claims (19)

1.圆柱阵列波浪力幅值波动特性、步长和包络线的获取方法，其特征在于，包括如下步骤：

步骤100，将波浪力幅值随无量纲波数变化形成的波浪力曲线中的多个高耸尖峰所在区域作为区域I，将高耸尖峰附近低于所述高耸尖峰且曲线波动间距随无量纲波数改变而变化的次级峰和谷所在区域作为区域II，将排除区域I和区域II后的波浪力曲线作为区域III；

步骤200，建立由排成一条直线的有限个相同坐底圆柱所组成的圆柱阵列坐标系并确定相关参数，将无量纲波数转化为相邻两圆柱间距离和波长的比值，则各圆柱的绕射波发生相长干涉时波浪力曲线出现峰点，发生相消干涉时波浪力曲线出现谷点；对于波浪入射角等于零的情形，分析入射波传到任意某个圆柱发生绕射及入射波传到该圆柱下游最后一个圆柱发生绕射这两个路径的波程差，求出波浪力曲线区域III中任意峰点横坐标的初步表达式和相邻峰点横坐标之差的初步表达式，以及任意谷点横坐标的初步表达式和相邻谷点横坐标之差的初步表达式，进而可知相邻两个峰点或相邻两个谷点的横坐标之差相等，最后得到任意圆柱波浪力曲线在区域III中的初步波动间距表达式；

步骤300，在任意某个圆柱上游处，将该圆柱产生的向上游传播的左传绕射波和该圆柱下游各个圆柱产生的左传绕射波求和，利用汉克尔函数取渐近表达式进行简化，得到对波动间距问题来说与当前圆柱阵列等价、且位置与现有圆柱位置均不相同的两个位于该圆柱下游的等价圆柱的具体位置，将离该圆柱更远的那个等价圆柱的位置代入初步波动间距表达式、任意峰点横坐标的初步表达式和任意谷点横坐标的初步表达式进行修正，即可得到波浪入射角等于零时波浪力曲线区域III中波动间距的最终表达式、任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式；

步骤400，对于波浪入射角不等于零的情形，利用前述同样方法，首先分析入射波传到任意某个圆柱发生绕射及入射波传到该圆柱上游的阵列端部第一个圆柱发生绕射这两个路径的波程差，然后分析入射波传到任意某个圆柱发生绕射及入射波传到该圆柱下游的最后一个圆柱发生绕射这两个路径的波程差，即可分别得到任意圆柱波浪力曲线在区域III中由于该圆柱上游圆柱绕射波的作用导致的任意峰点横坐标的初步表达式一、任意谷点横坐标的初步表达式一和初步波动间距表达式一，以及由于该圆柱下游圆柱绕射波的作用导致的任意峰点横坐标的初步表达式二、任意谷点横坐标的初步表达式二和初步波动间距表达式二；

步骤500，首先在任意某个圆柱的下游处，对该圆柱产生的向下游传播的右传绕射波和该圆柱上游各个圆柱产生的右传绕射波求和，利用前述同样的方法修正，得到波浪入射角不等于零时，任意峰点横坐标的最终表达式一、任意谷点横坐标的最终表达式一和最终波动间距表达式一；然后，在任意某个圆柱的上游处，对该圆柱产生的向上游传播的左传绕射波和该圆柱下游各个圆柱产生的左传绕射波求和，同样方法修正后得到任意峰点横坐标的最终表达式二、任意谷点横坐标的最终表达式二和最终波动间距表达式二；

步骤600，采用最终波动间距表达式一和最终波动间距表达式二给出的波动间距的更小者，就得到波浪入射角不等于零时任意圆柱波浪力曲线在区域III中最小波动间距的最终表达式；综合波浪入射角等于零时和不等于零时的波动间距的最终表达式，即可得到任意圆柱波浪力曲线在区域III中波动间距描述模型的最终表达式；

步骤700，根据入射角等于零和不等于零时波动间距的最终表达式计算出区域III中的最小波动间距，以该最小波动间距作为波浪力曲线在区域III中计算步长的上限，根据不同精度要求将最小波动间距除以相应的自然数作为区域III中计算步长的下限，从而得到区域III的计算步长；

根据所需精度要求将区域III的最小波动间距除以相应精度要求对应的自然数，得到区域II中的计算步长；根据所需精度要求将区域III的最小波动间距除以相应精度要求对应的自然数，得到区域I中的计算步长；

步骤800，采用最终表达式一和最终表达式二中横坐标的更小者作为波浪入射角不等于零时任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式，综合波浪入射角等于零时和不等于零时的最终表达式，即可得到任意圆柱波浪力曲线在区域III中任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式；根据横坐标的最终表达式得到任意峰点和谷点处的横坐标，从而可得对应的波数，求解线性方程组以获得速度势表达式中的未知绕射系数，进而可得到任意圆柱所受波浪力，将该波浪力无量纲化并取模，就得到波浪力曲线在区域III中任意峰点和谷点处的纵坐标；

根据上述获取的任意峰点和谷点处的横坐标以及任意峰点和谷点处的纵坐标，确定出每一个峰点和谷点的位置，将这些峰点连线，就得到波浪力曲线在区域III中的上包络线；将这些谷点连线，就得到波浪力曲线在区域III中的下包络线。

2.根据权利要求1所述的获取方法，其特征在于，

在所述圆柱阵列坐标系中，平面入射波传播方向与圆柱阵列整体坐标系中x轴正方向的夹角称为波浪入射角，且整体坐标系的建立使波浪入射角小于等于90度；k为圆柱阵列中的任意一根圆柱的编号，编号k增大方向与圆柱阵列整体坐标系中x轴正方向保持一致。

3.根据权利要求1所述的获取方法，其特征在于，

所述步骤200中的初步波动间距表达式、任意峰点横坐标的初步表达式和任意谷点横坐标的初步表达式的获取过程如下：

无量纲波数Kd/π可以改写为：Kd/π＝R/λ；

波浪力曲线中的第一个峰点和谷点对应的柱间距用下式表示：

2(N-k)R_p(1)＝λ

其中，K是波数，R＝2d是相邻圆柱的柱间距，λ是波长，脚标中，p代表峰点，v代表谷点，(1)代表第一个峰或者谷；

令在区域III中任意第s和第s+1个峰点对应的柱间距是R_p(s)和R_p(s+1)，R_p(s+1)＝R_p(s)+δR_p，峰点出现的条件是发生相长干涉，则相邻的第s和第s+1个峰点对应的波程差应该分别是s倍和s+1倍的波长，由此可得相邻峰点横坐标之差的初步表达式为：

谷点出现的条件是发生相消干涉，对于第s和第s+1个谷点对应的柱间距R_v(s)和R_v(s+1)＝R_v(s)+δR_v，相应的波程差分别等于2s-1倍和2s+1倍半波长，由此可得相邻谷点横坐标之差的初步表达式为：

由上述推导可知，相邻峰点和相邻谷点的横坐标之差相等，则初步波动间距表达式为：

N为圆柱阵列中的圆柱总数；此时波浪力曲线区域III中任意峰点横坐标的初步表达式和任意谷点横坐标的初步表达式为：

4.根据权利要求3所述的获取方法，其特征在于，

获取波浪入射角等于零时波浪力曲线区域III中波动间距的最终表达式、任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式的过程如下：

在任意k柱上游|x_k|处，其中x_k<0，把k柱以及k柱下游各柱的绕射势求和得到：

其中，

A_n为第1柱上的系数，i是虚数单位，ω是波浪圆频率，t是时间，n是整数，Z_n＝J′_n(Ka)/H′_n(Ka)，K为波数，J_n是第一类贝塞尔函数，H_n是第一类汉克尔函数，a是圆柱半径；

将汉克尔函数取渐近表达式进行简化后有：

其中，

上式第二项和第三项等价于两个与k柱左传波波程差为2(R/2)和2(N+1/2-k)R的等价圆柱左传波的叠加效果，这两个等价圆柱位于k柱下游R/2和(N+1/2-k)R处；

利用离k柱更远的等价圆柱的位置对初步波动间距表达式进行修正后即可得到波浪力曲线区域III中波动间距的最终表达式：

同样的，修正后的波浪力曲线区域III中任意峰点横坐标的最终表达式和任意谷点横坐标的最终表达式为：

5.根据权利要求3所述的获取方法，其特征在于，

对于波浪入射角不等于零的情况，任意峰点横坐标的最终表达式一、任意谷点横坐标的最终表达式一和最终波动间距表达式一的获取过程如下：

其中，k柱为陈列中的任意一个圆柱，β为波浪入射角，令β≠0时，在区域III中任意第s和第s+1个峰点、谷点对应的柱间距是

将这些表示柱间距的符号再增加上角标u后，表示由上游各柱右转波作用所导致的相应的量；

首先得到任意峰点横坐标的初步达式一、任意谷点横坐标的初步表达式一和初步波动间距表达式一：

然后对上述三式进行修正即得到任意峰点横坐标的最终表达式一、任意谷点横坐标的最终表达式一和最终波动间距表达式一；

其中“us”表示第k柱上游柱子的贡献。

6.根据权利要求5所述的获取方法，其特征在于，

对于波浪入射角不等于零的情况，任意峰点横坐标的最终表达式二、任意谷点横坐标的最终表达式二和最终波动间距表达式二的获取过程如下：

其中，

等表示柱间距的符号再增加上角标l表示由下游各柱左传波作用所导致的相应的量；

首先得到任意峰点横坐标的初步表达式二、任意谷点横坐标的初步表达式二和初步波动间距表达式二：

然后对上述三式进行修正，即得到任意峰点横坐标的最终表达式二、任意谷点横坐标的最终表达式二和最终波动间距表达式二：

其中下标“ds”表示第k柱下游柱子的贡献。

7.根据权利要求6所述的获取方法，其特征在于，

所述步骤600中的最终表达式在波浪入射角等于零时，任意圆柱波浪力曲线在区域III中波动间距为：

在波浪入射角不等于零时，任意圆柱波浪力曲线在区域III中波动间距为：

8.根据权利要求1所述的获取方法，其特征在于，

所述区域III的自然数取值范围在2～10之间，所述区域II的自然数取值范围在5～10之间，所述区域I的自然数取值范围在40～50之间；当所述区域III的计算步长取五分之一的区域III最小波动间距，所述区域II的计算步长取十分之一的区域III最小波动间距，所述区域I的计算步长取五十分之一的区域III最小波动间距时，波浪力曲线的计算精度达到相对误差在1％以内。

9.根据权利要求6所述的获取方法，其特征在于，

获取各所述峰点的横坐标的过程如下：

对于波浪入射角等于零的情况，波浪力曲线区域III中任意峰点横坐标的以柱间距-波长比表示的最终表达式为:

对于波浪入射角不等于零的情况，结合任意峰点横坐标的最终表达式一和任意峰点横坐标的最终表达式二可以得到以柱间距-波长比表示的任意峰点横坐标的最终表达式：

获取各所述谷点的横坐标的过程如下：

对于波浪入射角等于零的情况，波浪力曲线区域III中任意谷点横坐标的以柱间距-波长比表示的最终表达式为:

对于波浪入射角不等于零的情况，结合任意谷点横坐标的最终表达式一和任意谷点横坐标的最终表达式二可得到以柱间距-波长比表示的任意谷点横坐标的最终表达式：

10.根据权利要求9所述的获取方法，其特征在于，

每一个所述峰点坐标和所述谷点坐标中的纵坐标获取方式如下：

根据坐底圆柱阵列水波绕射问题中任意k柱附近速度势的空间因子φ(r_k,θ_k)公式：

其中，

为绕射系数，k为圆柱阵列中任意一根圆柱的编号，编号k增大方向与x轴正方向保持一致，(r_k,θ_k)为垂直轴z轴通过k柱轴线的局部圆柱坐标系的极坐标，Z_n＝J′_n(Ka)/H′_n(Ka)，K为波数，a为圆柱半径，J_n为第一类贝塞尔函数，H_n是第一类汉克尔函数，n为整数；

求解所述速度势表达式中绕射系数的线性方程组如下：

其中，β为平面入射波传播方向与圆柱阵列整体坐标系中x轴正方向的夹角，即为波浪入射角，且整体坐标系的建立使波浪入射角β≤π/2，R_jk为第k个柱子到第j个柱子的距离，i为虚数单位，m为整数,α_jk为第k个柱子到第j个柱子的方向角，I_k为入射波在第k个柱子的相位因子；

将所述峰点和所述谷点对应的波数K代入该方程可求得所述峰点和所述谷点对应波数下绕射系数

的值，将绕射系数

代入下式

就可得到波浪力曲线在区域III中任意峰点和谷点对应波数下，任意圆柱k所受沿圆柱阵列水平截面中各圆心连线方向的波浪力F^k；其中，ρ为水的密度，g为重力加速度，A为入射波的波幅，h为水深；

用相同几何尺寸的圆柱在相同环境条件下所受波浪力对前式所示的波浪力进行无量纲化，可得到波浪力曲线在区域III中任意峰点和谷点对应波数下圆柱阵列中任意第k柱的无量纲波浪力为：

对该式取模即得到无量纲波浪力幅值，这也就是波浪力曲线在区域III中任意峰点和谷点处的纵坐标。

11.根据权利要求6所述的获取方法，其特征在于，

所述步骤200至步骤600中描述模型的获取过程还可采用如下获取方法：

步骤201，建立由排成一条直线的有限个相同坐底圆柱所组成的圆柱阵列的整体坐标系，列出坐底圆柱阵列水波绕射问题中任意圆柱附近的速度势表达式，速度势表达式中由线性方程组确定的未知系数称为绕射系数，由于绕射系数幅值随无量纲波数变化形成的绕射系数曲线与波浪力幅值随无量纲波数变化形成的波浪力曲线的波动特性和波动间距一致，因此对波浪力曲线的研究转化为对绕射系数曲线的研究；针对确定绕射系数的线性方程组采用克莱姆法则求解，用上述线性方程组右端列向量替换该线性方程组左端系数方阵相应的列然后得到目标方阵行列式；

步骤301，忽略高阶小量得到目标方阵行列式模的平方表达式；忽略平方表达式中随无量纲波数缓慢变化的慢变函数，得到与波动间距分析密切相关的目标表达式；对目标表达式求一阶导数，令该一阶导数为零，得到超越方程；

步骤401，利用超越方程分别推导波浪入射角等于零和不等于零时任意圆柱波浪力曲线在区域III中的波动间距表达式，综合各表达式后即得到区域III的描述模型。

12.根据权利要求11所述的获取方法，其特征在于，

所述步骤201中的速度势表达式如下：

其中，

为绕射系数，k为圆柱阵列中任意一根圆柱的编号，编号k增大方向与x轴正方向保持一致，(r_k,θ_k)为垂直轴z轴通过k柱轴线的局部圆柱坐标系的极坐标，Z_n＝J′_n(Ka)/H′_n(Ka)，K为波数，a为圆柱半径，J_n为第一类贝塞尔函数，H_n是第一类汉克尔函数，n为整数，i为虚数单位。

13.根据权利要求12所述的获取方法，其特征在于，

求解所述速度势表达式中绕射系数的线性方程组如下：

其中，β为波浪入射角，是由平面入射波传播方向与圆柱阵列整体坐标系中x轴正方向形成的夹角，且整体坐标系的建立使波浪入射角β≤π/2，K为波数,R_jk为第k个柱子轴线到第j个柱子轴线的距离，i为虚数单位，m为整数,N为圆柱阵列中的圆柱总数，α_jk为第k个柱子到第j个柱子的方向角，I_k为入射波在第k个柱子的相位因子；

利用克莱姆法则求解上述线性方程组后得到的绕射系数如下：

其中，D为所述线性方程组中系数方阵的行列式，

为用所述线性方程组的右端列向量替换系数方阵对应列得到的目标方阵行列式。

14.根据权利要求13所述的获取方法，其特征在于，

对绕射系数

的分析可进一步转化为对目标方阵行列式模的平方

的分析，将

展开并忽略高阶小量，进而求其模的平方后得到的公式即为目标方阵行列式模的平方表达式：

其中，

κ＝Kd/π为无量纲波数，2d为相邻圆柱轴线之间的距离，M为截断项数。

15.根据权利要求14所述的获取方法，其特征在于，

所述目标表达式为：

对目标表达式求一阶导数，令该一阶导数为零，得到的所述超越方程为：

[2(N-k)+1]tan(2πκ)＝tan{[2(N-k)+1]2πκ}。

16.根据权利要求15所述的获取方法，其特征在于，

利用超越方程推导任意圆柱k在波浪入射角等于零时其在区域III中的波浪力曲线波动间距的表达式过程为：

步骤410，对超越方程左右两端表达式分别做泰勒展开以得到各自的近似表达式；

步骤411，再把各近似表达式分别代入超越方程，解析求得任意圆柱k在区域III波浪力曲线波动间距中的上下限表达式；

步骤412，忽略小量后，上述上下限表达式相同，因而可得到波浪入射角等于零时，圆柱阵列中任意圆柱k波浪力曲线在区域III中的波动间距描述模型的表达式为

17.根据权利要求15所述的获取方法，其特征在于，

利用超越方程推导任意圆柱k在波浪入射角不等于零时波浪力曲线在区域III中的波动间距表达式时，需要将超越方程[2(N-k)+1]tan(2πκ)＝tan{[2(N-k)+1]2πκ}中的2πκ用(1+cosβ)πκ、(1-cosβ)πκ代替以得到波浪入射角不等于零时的两个修正超越方程，然后对这两个修正超越方程的左、右两端表达式分别做泰勒展开得到各自近似表达式；再把近似表达式分别代入两个修正超越方程，则可解析求得区域III波浪力曲线波动间距上下限的表达式，进而可以得到波浪入射角不等于零情况下的波浪力曲线波动间距表达式。

18.根据权利要求17所述的获取方法，其特征在于，

在波浪入射角不等于零时，波浪力曲线按波动特性以及波动间距不同分为以下五种情形：

(1)[2(N-k)+1](1+cosβ)＞＞[2(k-1)+1](1-cosβ)；

(2)[2(k-1)+1](1-cosβ)＞＞[2(N-k)+1](1+cosβ)；

(3)[2(N-k)+1](1+cosβ)>[2(k-1)+1](1-cosβ)＞＞4；

(4)[2(k-1)+1](1-cosβ)>[2(N-k)+1](1+cosβ)＞＞4；

(5)[2(N-k)+1](1+cosβ)和[2(k-1)+1](1-cosβ)量级比较接近，均远大于4；

其中的“＞＞4”和“远大于4”表示区域III的波动间距远小于区域I各个高耸尖峰之间的距离，

第(1)种情形的波浪力曲线区域III波动间距表达式为：

第(2)种情形的波浪力曲线区域III波动间距表达式为：

第(3)种情形的波浪力曲线区域III波动间距表达式采用第(1)种情形和第(2)种情形间距表达式中的波动间距更小者；

第(4)种情形的波浪力曲线区域III波动间距表达式采用第(1)种情形和第(2)种情形间距表达式中的波动间距更小者；

第(5)种情形的波浪力曲线区域III波动间距表达式采用第(1)种情形和第(2)种情形间距表达式中的波动间距更小者。

19.根据权利要求18所述的获取方法，其特征在于，

综合在波浪入射角不等于零时五种情形下区域III的波动间距表达式，得到波浪入射角不等于零时圆柱阵列中任意圆柱k波浪力曲线在区域III中最小波动间距

描述模型的表达式为：

Images