CN109726419A - 基于相长相消的圆柱阵列波浪力曲线计算步长的确定方法 - Google Patents

Abstract

本发明提供了基于相长相消的圆柱阵列波浪力曲线计算步长的确定方法，通过建立由排成一条直线的有限个相同坐底圆柱所组成的圆柱阵列坐标系并确定相关参数，分析波浪入射角等于零和不等于零时各圆柱的左转绕射波和右转绕射波情况，得到任意圆柱波浪力曲线在非捕获区中波动间距描述模型的最终表达式；再以此为基础确定各区域的计算步长。本发明的预测公式可以在保证精度的前提下减少工作量；本发明的计算步长确定方法可以在不盲目增加不必要计算时间的前提下，得到准确的波浪力曲线。

Description

基于相长相消的圆柱阵列波浪力曲线计算步长的确定方法

技术领域

本发明涉及海洋工程领域，特别涉及由大数量穿透水面圆柱组成的圆柱阵列在波浪作用下，利用波的相长相消原理建立任意圆柱所受波浪力的波浪力幅值曲线计算步长的确定方法。

背景技术

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

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

如图1所示的单排坐底圆柱阵列中任意圆柱上所受波浪力的幅值是随无量纲波数上下起伏、波动变化的。一般来说，单排较大数量(例如，数量大于9)有限圆柱阵列中单个柱子所受到波浪力的幅值随波数变化具有三个明显的特点：1)波浪力幅值随无量纲波数变化形成的波浪力曲线具有几个高耸的尖峰，这几个高耸尖峰所在区域本发明称为区域I(Region I)；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，建立由排成一条直线的有限个相同坐底圆柱所组成的圆柱阵列坐标系并确定相关参数，将无量纲波数转化为相邻两圆柱间距离和波长的比值；则确定波浪力曲线峰点和谷点的方法可描述为：当各圆柱的绕射波发生相长干涉时波浪力曲线出现峰点，当各圆柱绕射波发生相消干涉时波浪力曲线出现谷点；对于波浪入射角等于零的情形，分析入射波传到任意某个圆柱发生绕射及入射波传到该圆柱下游最后一个圆柱发生绕射这两个路径的波程差，该波程差与相邻两圆柱间距离成正比，令该波程差分别等于s倍和s+1倍的波长，可求得波浪力曲线区域III中任意峰点的横坐标以及相邻峰点横坐标之差的初步表达式，令该波程差分别等于2s-1倍和2s+1倍半波长，可求得波浪力曲线区域III中任意谷点的横坐标以及相邻谷点横坐标之差的初步表达式；进而可知相邻两个峰点或相邻两个谷点的横坐标之差相等，因此可得到任意圆柱波浪力曲线在区域III中的初步波动间距表达式；

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

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

步骤500，首先在任意某个圆柱的下游处，对该圆柱产生的向下游传播的右传绕射波和该圆柱上游各个圆柱产生的右传绕射波求和，利用与步骤300同样的方法修正，即可得到波浪入射角不等于零时，任意圆柱波浪力曲线在区域III中由于该圆柱上游圆柱绕射波的作用导致的任意峰点和谷点的横坐标以及最终波动间距表达式一；然后，在任意某个圆柱的上游处，对该圆柱产生的向上游传播的左传绕射波和该圆柱下游各个圆柱产生的左传绕射波求和，利用步骤300同样的方法修正，即可得到波浪入射角不等于零时，任意圆柱波浪力曲线在区域III中由于该圆柱下游圆柱绕射波的作用导致的任意峰点和谷点的横坐标以及最终波动间距表达式二；

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

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

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

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

在本发明的一个实施方式中，所述步骤200中的初步表达式获取过程如下：

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

波浪力曲线中的第一个峰点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倍的波长，由此可得相邻峰点横坐标之差的初步表达式为：

考虑发生相消干涉的条件，对于第s和第s+1个谷点对应的柱间距R_v(s)和R_v(s+1)＝R_v(s)+δR_v，可得相邻谷点横坐标之差的初步表达式为：

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

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柱上游的圆柱阵列端部第一个圆柱向下游传播的右传绕射波与k柱右传绕射波的波程差，使其满足发生相长或相消干涉的条件，可得到初步波动间距表达式一，过程如下：

其中，k柱为圆柱阵列中的任意一个圆柱，β为波浪入射角，令β≠0时，在区域III中任意第s和第s+1个峰点、谷点对应的柱间距是 将这些表示柱间距的符号再增加上角标u后，表示由上游各柱右传波作用所导致的相应的量；

得到初步波动间距表达式一：

在本发明的一个实施方式中，对所述初步波动间距表达式一进行修改，得到波浪入射角不等于零时最终波动间距表达式一及此状态下的任意峰点或谷点的横坐标的最终表达式一分别为；

在本发明的一个实施方式中，对于波浪入射角不等于零的情况，再分析位于k柱下游的最后一个圆柱向上游传播的左传绕射波与k柱左传绕射波的波程差，使其满足发生相长或相消干涉的条件，可得到初步波动间距表达式二，过程如下：

其中，等表示柱间距的符号，再增加上角标l表示由下游各柱左传波作用所导致的相应的量；得到初步波动间距表达式二如下：

在本发明的一个实施方式中，对所述初步波动间距表达式二进行修改，得到波浪入射角不等于零时最终波动间距表达式二及此状态下任意峰点或谷点的横坐标的最终表达式二分别为：

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

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

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

本发明提供了这样的认识和理解：波浪力幅值随无量纲波数变化形成的波浪力曲线中非捕获区(区域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的波动间距测量值和理论计算值随柱子编号变化的图像。

具体实施方式

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

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

目前已有不少文献研究了无限长圆柱阵列或水槽中心线上布置单个圆柱的trapped mode频率问题，这些结果可以估算有限数目坐底圆柱阵列的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，建立由排成一条直线的有限个坐底圆柱所组成的圆柱阵列坐标系并确定相关参数，将无量纲波数转化为相邻两圆柱间距离和波长的比值；则确定波浪力曲线峰点和谷点的方法可描述为：当各圆柱的绕射波发生相长干涉时波浪力曲线出现峰点，当各圆柱绕射波发生相消干涉时波浪力曲线出现谷点；对于波浪入射角等于零的情形，分析入射波传到任意某个圆柱发生绕射及入射波传到该圆柱下游最后一个圆柱发生绕射这两个路径的波程差，该波程差与相邻两圆柱间距离成正比，令该波程差分别等于s倍和s+1倍的波长(s为任意自然数)，可求得波浪力曲线区域III中任意峰点的横坐标以及相邻峰点横坐标之差的初步表达式，令该波程差分别等于2s-1倍和2s+1倍半波长，可求得波浪力曲线区域III中任意谷点的横坐标以及相邻谷点横坐标之差的初步表达式；进而可知相邻两个峰点或相邻两个谷点的横坐标之差相等，因此可得到任意圆柱波浪力曲线在区域III中的初步波动间距表达式；

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

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

如图6所示，其中的初步表达式获取过程如下：

无量纲波数(即波浪力曲线横坐标)Kd/π可以改写为：Kd/π＝R/λ；于是，波浪力幅值随无量纲波数变化也可以认为是波浪力幅值随着柱间距—波长比而变化。当j柱的绕射波传到k柱附近与k柱绕射波发生相长(相消)干涉时，k柱上的波浪力幅值会得到一个峰(谷)值。如图6所示，R_p(1)和R_v(1)分别表示波浪力曲线中的第一个峰点和谷点对应的柱间距，对于波浪入射角等于零的情况，波浪力曲线中的第一个峰点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，可以有同样的推导：

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

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

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

步骤300，在任意某个圆柱上游处，将该圆柱产生的向上游传播的左传绕射波和该圆柱下游各个圆柱产生的左传绕射波求和，利用汉克尔函数取渐近表达式进行简化，得到对波动间距问题来说与当前圆柱阵列等价、且位置与现有圆柱位置均不相同的两个位于该圆柱下游的等价圆柱的具体位置，将离该圆柱更远的那个等价圆柱的位置代入初步波动间距表达式对该表达式进行修正，即可得到波浪入射角等于零时波浪力曲线区域III中波动间距的最终表达式、和此时波浪力曲线区域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，对于波浪入射角不等于零的情形，利用步骤200同样的方法，首先分析入射波传到任意某个圆柱发生绕射及入射波传到该圆柱上游的阵列端部第一个圆柱发生绕射这两个路径的波程差，然后分析入射波传到任意某个圆柱发生绕射及入射波传到该圆柱下游的最后一个圆柱发生绕射这两个路径的波程差，即可分别得到任意圆柱波浪力曲线在区域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，首先在任意某个圆柱的下游处，对该圆柱产生的向下游传播的右传绕射波和该圆柱上游各个圆柱产生的右传绕射波求和，利用步骤300同样的方法修正，即可得到波浪入射角不等于零时，任意圆柱波浪力曲线在区域III中由于该圆柱上游圆柱绕射波的作用导致的任意峰点和谷点的横坐标以及最终波动间距表达式一；然后，在任意某个圆柱的上游处，对该圆柱产生的向上游传播的左传绕射波和该圆柱下游各个圆柱产生的左传绕射波求和，利用步骤300同样的方法修正，即可得到波浪入射角不等于零时，任意圆柱波浪力曲线在区域III中由于该圆柱下游圆柱绕射波的作用导致的任意峰点和谷点的横坐标以及最终波动间距表达式二；

与β＝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时，单排坐底圆柱群阵列第k柱受到的波浪力在区域III的振荡间距测量值和描述模型表达式计算值的对比。通过对比发现，本发明的描述模型表达式预测值与实际计算值吻合非常好。

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

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

在入射角等于零时(β＝0)的最终表达式

和

入射角不等于零时(β≠0)的最终表达式

中代入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％以内。

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

Claims (10)

1.基于相长相消的圆柱阵列波浪力曲线计算步长的确定方法，其特征在于，包括如下步骤：

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

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

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

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

步骤500，首先在任意某个圆柱的下游处，对该圆柱产生的向下游传播的右传绕射波和该圆柱上游各个圆柱产生的右传绕射波求和，利用与步骤300同样的方法修正，即可得到波浪入射角不等于零时，任意圆柱波浪力曲线在区域III中由于该圆柱上游圆柱绕射波的作用导致的任意峰点和谷点的横坐标以及最终波动间距表达式一；然后，在任意某个圆柱的上游处，对该圆柱产生的向上游传播的左传绕射波和该圆柱下游各个圆柱产生的左传绕射波求和，利用步骤300同样的方法修正，即可得到波浪入射角不等于零时，任意圆柱波浪力曲线在区域III中由于该圆柱下游圆柱绕射波的作用导致的任意峰点和谷点的横坐标以及最终波动间距表达式二；

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

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

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

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

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

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

所述步骤200中的初步表达式获取过程如下：

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

波浪力曲线中的第一个峰点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倍的波长，由此可得相邻峰点横坐标之差的初步表达式为：

考虑发生相消干涉的条件，对于第s和第s+1个谷点对应的柱间距R_v(s)和R_v(s+1)＝R_v(s)+δR_v，可得相邻谷点横坐标之差的初步表达式为：

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

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.根据权利要求1所述的确定方法，其特征在于，

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

其中，k柱为圆柱阵列中的任意一个圆柱，β为波浪入射角，令β≠0时，在区域III中任意第s和第s+1个峰点、谷点对应的柱间距是 将这些表示柱间距的符号再增加上角标u后，表示由上游各柱右传波作用所导致的相应的量；

得到初步波动间距表达式一：

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

对所述初步波动间距表达式一进行修改，得到波浪入射角不等于零时最终波动间距表达式一及此状态下的任意峰点或谷点的横坐标的最终表达式一分别为:

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

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

其中， 等表示柱间距的符号，再增加上角标l表示由下游各柱左传波作用所导致的相应的量；得到初步波动间距表达式二如下：

8.根据权利要求7所述的确定方法，其特征在于，

对所述初步波动间距表达式二进行修改，得到波浪入射角不等于零时最终波动间距表达式二及此状态下任意峰点或谷点的横坐标的最终表达式二分别为：

9.根据权利要求8所述的确定方法，其特征在于，

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

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

10.根据权利要求1所述的确定方法，其特征在于，

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