CN109726422A

CN109726422A - 基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法

Info

Publication number: CN109726422A
Application number: CN201810784158.0A
Authority: CN
Inventors: 曾晓辉; 于法军; 孙哲; 史民
Original assignee: Institute of Mechanics of CAS; University of Chinese Academy of Sciences
Current assignee: Institute of Mechanics of CAS; University of Chinese Academy of Sciences
Priority date: 2018-07-17
Filing date: 2018-07-17
Publication date: 2019-05-07
Anticipated expiration: 2038-07-17
Also published as: CN109726422B

Abstract

本发明提供了一种基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法，包括如下步骤：列出任意圆柱附近水波绕射的速度势，其中绕射系数由线性方程组确定；采用克莱姆法则求解该线性方程组得到目标方阵行列式；忽略该行列式模的平方表达式中随波数缓慢变化的量，得到目标表达式；对该表达式求一阶导数并令其等于零，得到超越方程后推导出非捕获区任意圆柱波浪力幅值曲线的波动间距表达式，根据所需精度除以相应的自然数，即可得到不同区域的计算步长。本发明阐明波浪力幅值曲线在非捕获区中波动间距不随无量纲波数的变化而改变，只与圆柱阵列中圆柱总数、标识圆柱位置的柱子编号和波浪入射角有关，同时可以在保证精度的前提下减少计算工作量。

Description

基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法

技术领域

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

背景技术

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

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

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

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

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

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

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

发明内容

本发明的目的是要提供一种由较大数量穿透水面圆柱组成的圆柱阵列在波浪作用下，基于绕射理论的任意圆柱所受波浪力的波浪力幅值曲线计算步长的确定方法。

特别地，本发明提供一种基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法，包括如下步骤：

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

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

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

步骤400，利用超越方程分别推导波浪入射角等于零和不等于零时任意圆柱波浪力曲线在区域III中的波动间距表达式，然后根据该波动间距表达式计算出区域III中的最小波动间距，以该最小波动间距作为波浪力曲线在区域III中计算步长的上限，根据不同精度要求将最小波动间距除以相应的自然数作为区域III中计算步长的下限，从而得到区域III的计算步长；

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

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

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

其中，为绕射系数，整体坐标系x轴穿过圆柱阵列水平截面中的圆心连线，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)波动特性的理解和认识，本发明提供了非捕获区中波动间距的预测公式。基于本发明，在进行有关工程结构设计和评估时，可以在保证精度的前提下减少工作量，缩短设计和评估周期，为提高工程结构的设计和评估水平提供技术支撑。

采用本发明给出的波浪力曲线计算步长确定方法，只要已知圆柱阵列中圆柱总数、标识圆柱位置的圆柱编号和波浪入射角，就可以在任意圆柱波浪力曲线计算开始之前，准确预测出计算步长取值多大即可得到准确的波浪力曲线。从而可以在不盲目增加不必要计算时间的前提下，得到准确的波浪力曲线。

附图说明

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

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

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

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

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

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

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

图10为圆柱总数为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的计算范围。对于所得到尖峰对应的波数值，绘制出波浪力幅值与柱子编号的关系图，如能呈现一个完整的半波形式，如图4所示，最大的波浪力作用在中间柱上，则说明该尖峰位置足够准确。如果不是如此，还可以继续加密步长找到更准确的尖峰点。

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

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

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

圆柱数目为N/3的单排坐底圆柱阵列之中间柱波浪力曲线尖峰位置横坐标对应着圆柱数目为N的圆柱阵列之中间柱波浪力曲线尖峰左边第一个峰点位置的横坐标，此峰点位置横坐标对应的无量纲波数下圆柱数目为N的圆柱阵列的波浪力幅值与圆柱编号关系图呈现出三个半波的形式，如图6所示，这三个半波的最高峰对应的波浪力幅值与同样波数下圆柱数目为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柱的情况，可以参照上述中间柱的范围来确定。

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

如图1所示，本发明一个实施例的基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法，包括如下步骤：

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

步骤200，建立由排成一条直线的有限个相同坐底圆柱所组成的圆柱阵列的整体坐标系，列出坐底圆柱阵列水波绕射问题中任意圆柱附近的速度势表达式，速度势表达式中由线性方程组确定的未知系数称为绕射系数，由于绕射系数幅值随无量纲波数变化形成的绕射系数曲线与波浪力幅值随无量纲波数变化形成的波浪力曲线的波动特性和波动间距一致，因此对波浪力曲线的研究转化为对绕射系数曲线的研究；针对确定绕射系数的线性方程组采用克莱姆法则(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)，a为圆柱半径，J_n为第一类贝塞尔函数，H_n是第一类汉克尔函数，n为整数；β为波浪入射角，是由平面入射波传播方向与圆柱阵列整体坐标系中x轴正方向形成的夹角，且整体坐标系的建立使波浪入射角β≤π/2，K为波数,R_jk为第k个柱子轴线到第j个柱子轴线的距离，i为虚数单位，m为整数,N为圆柱阵列中的圆柱总数，α_jk为第k个柱子到第j个柱子的方向角，I_k为入射波在第k个柱子的相位因子。

开展有关计算和分析后，可以发现波浪力曲线区域III与绕射系数幅值随无量纲波数变化形成的绕射系数曲线的区域III范围相同，波动间距也相同。因此对波浪力曲线波动间距的讨论可以转化为对绕射系数随无量纲波数变化形成的绕射系数曲线的波动间距的研究。为了得到波动间距的描述模型，我们依据克莱姆法则(Cramer's Rule)求解绕射系数线性方程组，绕射系数可表示为

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

经过研究发现，对绕射系数的分析可进一步转化为对目标方阵行列式模的平方的分析，将展开并忽略高阶小量，进而求其模的平方后再忽略高阶小量，得到的公式即为目标方阵行列式模的平方表达式：

其中，

上述式中，κ＝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)

以下说明波浪入射角等于零(β＝0)和不等于零(β≠0)两种情况下的区域III波动间距的描述模型和表达式。

一、波浪入射角等于零时，单排坐底圆柱阵列波浪力曲线区域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…附近。图7(a)和(b)给出了交点的分布情况。显然，交点均在χ(κ)的渐近线附近。为了得到交点横坐标的解析表达式，首先把χ(κ)在qπ+π/2附近展开。于是在qπ+π/2，q＝0,1,2…附近有

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

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

如图7(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个解

如图7(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)，波浪力曲线最小波动间距描述模型的表达式为：

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

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

通过区域III的波动间距表达式可计算出区域III中的最小波动间距，即对

代入N，k和β后即可计算出最小波动间距或以下为简洁起见，统一用符号表示波浪入射角等于零和不等于零两种情况下的最小波动间距，需注意的是波浪入射角β＝0时采用式(23)计算最小波动间距；此时的最小波动间距作为区域III计算步长的上限，而计算步长的下限，则根据计算时的精度要求选取一个范围在2～10之间的自然数，作为最小波动间距的除数，即可确定该下限，这里自然数的取值越大，则精度越高，相应花费的计算时间也越长。通过确定后的上限和下限，即可得到区域III的计算步长。

由于区域II范围比区域III范围小很多且波动间距变小，因此，可根据精度要求取一个范围在5～10之间的自然数作为区域III的最小波动间距的除数，即可得到区域II的计算步长，这里同样是自然数的取值越大，则精度越高，相应花费的计算时间也越长。

由于区域I范围比区域II小很多且具有分离的高耸尖峰，因此，可根据精度要求取一个范围在40～50之间的自然数作为区域III的最小波动间距的除数，即可得到区域I的计算步长，这里同样是自然数的取值越大，则精度越高，相应花费的计算时间也越长。

具体地，当区域III的计算步长取五分之一的区域III最小波动间距区域II的计算步长取十分之一的区域III最小波动间距区域I的计算步长取五十分之一的区域III最小波动间距时，波浪力曲线的计算精度达到相对误差在1％以内。

图9是柱子总数N＝101,波浪入射角β＝0，直径柱间距比a/d＝1/4时，波浪力幅值的波动间距测量值随着无量纲波数变化的曲线。可以看到，波动间距在很大波数范围下是恒定的，这个区域就是区域III，之后波动间距迅速下降的区域就是区域I和II，图中渐近线及数值是利用本发明最终表达式计算得到的理论预测值，结果非常吻合。

图10分别是N＝11,N＝21,N＝51,N＝101情况下，波浪入射角β＝0时，单排坐底圆柱群阵列第k柱受到的波浪力在区域III的波动间距测量值和描述模型表达式计算值的对比。通过对比发现，本发明的描述模型表达式预测值与实际计算值吻合非常好。

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

Claims

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

2.根据权利要求1所述的基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法，其特征在于，

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

3.根据权利要求1所述的基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法，其特征在于，

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

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

4.根据权利要求3所述的基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法，其特征在于，

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

5.根据权利要求4所述的基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法，其特征在于，

对绕射系数的分析可进一步转化为对目标方阵行列式模的平方的分析，将展开并忽略高阶小量，进而求其模的平方后得到的公式即为目标方阵行列式模的平方表达式：

其中，

6.根据权利要求5所述的基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法，其特征在于，

所述目标表达式为：

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

7.根据权利要求6所述的基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法，其特征在于，

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

8.根据权利要求6所述的基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法，其特征在于，

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

9.根据权利要求8所述的基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法，其特征在于，

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

(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；

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

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

10.根据权利要求9所述的基于绕射理论的圆柱阵列波浪力曲线计算步长确定方法，其特征在于，

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