多个椭圆柱波浪力的一种解析解1)

赵 密 龙彭振 王丕光，2) 张 超 杜修力

* (北京工业大学城建学部建筑工程学院，北京 100124)

† (福州大学土木工程学院，福州 350116)

引言

近年来，中国的近海工程迅猛增加，如人工岛，跨海桥梁等，对近海结构的研究日益受到人们的重视.由于上述工程常年遭受波浪的威胁，其支撑结构可能承受着强度相当大的波浪力，而支撑结构广泛使用了多柱体，柱体截面形式不再局限于圆形，因此对非圆多柱体体系中波浪力的研究可为实际工程提供一定的理论指导.

近海工程中对波浪力的研究根据结构尺寸，可以划分为对小尺度结构所受波浪力的研究和对大尺度结构所受波浪力的研究.其中小尺度结构上的波浪力可以通过莫里森方程求解[1];当结构截面尺寸与波长之比大于0.2 时，大尺寸结构阻碍了波浪运动，波浪在结构表面生成不可忽略的散射波，因此波浪力的计算需要考虑波浪散射现象[2-3].MacCamy和Fuchs[4]提出了绕射波浪理论来计算有限水域中大尺寸圆柱线性波浪力.当波陡大时，波浪中的非线性成分增多，用线性波浪理论计算的结果误差增大，国内外学者就此提出了解决波浪二阶作用的方法[5-6].Tao 等[7]和Song 等[8]利用比例边界元法(SBFEM)分析短峰波在圆柱上绕射以及多个任意截面形状柱体和线性波浪相互作用，比例边界元对于解决波浪在柱面绕射问题较高的求解效率.Wang 等[9]利用有限元方法，解决了倾斜圆柱的波浪绕射问题.流场中只存在单个结构的情形很少，大多数情况下近海结构的承重构件由多柱体构成.Linton 和Evans[10]提出了计算多个圆柱附近波浪速度势的简化公式.Kagemoto 和Yue[11]提出了一种精确代数方法来分析水波中多个三维物体之间的相互作用问题.Chen 等[12]利用Null-Field 积分方程来解决波浪在多直立圆柱体中散射问题，以免带来用格林函数分析产生的不规则频率问题.缪国平等[13]假设柱体表面有若干波动源，就100 根以上的单排圆桩体系在波浪中的响应进行了研究，得到了多柱体对流体干扰的特点.Wang 等[14]将流场分为结构附近的有限域和无限域，用椭圆吸收边界模拟无限域，最后利用有限元方法求出波浪对群桩的影响，该方法可应用到任意光滑截面的多柱体上.

实际工程中同样存在非圆形截面柱体，如椭圆形和类椭圆等.Chen 和Mei[15]基于椭圆柱坐标系引入Mathieu 函数，将椭圆形固定平台表面的波浪力分解为角向分量和径向分量，分析了结构表面绕射问题，并得到相应解.William[16]简化了文献[15]中复杂级数展开，并得到两种近似方法，第一种适用于小曲率椭圆截面柱体，另一种利用格林函数获得了柱体表面的流体速度势.Liu 等[17-18]利用线性绕射理论和傅里叶展开，推导了任意光滑截面的铅直柱体和截断柱体表面波浪力.Zhang 和Williams[19]研究了水平刚性椭圆形薄板对入射波浪的绕射效应.Bhatta[20]推导了椭圆截面桩和圆形截面桩的波浪力计算公式，并实现了前者到后者的变换.Wang 等[21-22]首先提出针对单根椭圆柱上地震引起的动水力的解析解，再利用有限元方法改进该发法为适用于工程的附加质量矩阵方法，并且针对短峰波浪下的椭圆柱，推导得到波浪力的解析表达式，比较短峰波和一般线性波下椭圆柱周围的波爬.

对于多个非圆形截面柱体阵列结构，Chatjigeorgiou等[23-25]利用Bessel 函数和Mathieu 函数的加法定理，依次研究了流体中椭圆截面群桩波浪散射问题，椭圆柱和圆柱间的波浪散射问题和截断椭圆群桩上的波浪散射问题.Lee[26]则应用Collocation Multipole方法来解决椭圆截面多柱体间散射问题，免除了Mathieu 函数加法定理的繁复.Liu 等[27]利用傅里叶级数将径向函数展开，以求解任意光滑截面群桩(如余弦形截面、类椭圆截面等)上的波浪力.

上述多个柱体波浪力解析分析方法中，多数考虑了入射波及入射波引起的单柱体散射波，但没有考虑该散射波入射其他柱体后产生的第二次散射波以及高次散射波[25].本文基于椭圆坐标系和绕射波理论等，推导得到考虑高次散射波的椭圆多柱体波浪力公式，研究波数等参数变化时高次波对总波浪力的贡献，以期为实际工程中估算波浪力提供参考数据.

1 问题描述

1.1 控制方法和边界条件

流场中椭圆形截面柱体的阵列如图1 所示.假设所有柱体刚性，a为椭圆的半长轴，b为椭圆的半短轴，底端固定在深度为h的流体中，流体无旋、无黏性并且不可压缩.该多柱体体系受到与x轴成α，频率为 ω，幅值为H/2 的简谐波浪作用.

图1 柱体布置和坐标系统Fig.1 Arrangement of bodies and coordinates systems

本文使用的椭圆柱坐标系(ξ，μ，z)中，ξ和μ分别为椭圆柱坐标系的径向坐标和角向坐标，两坐标互相正交;坐标原点位于椭圆柱底部，z轴沿柱体轴线向上.椭圆坐标系和笛卡尔直角坐标系转换关系如下

式中，μ2=a2-b2；ξ，μ和z取值范围分别为0 ≤ξ＜∞，0≤μ＜2π和0≤z≤h.

简谐波浪下流场中波浪压力为

椭圆柱坐标系中，频域下的波浪压力p满足如下流体控制方程

水体自由表面、刚性地面、水体与结构交界面和无穷远处的边界条件分别为

式中，波数k满足 ω2=gktanh(kh) ; ξ0为柱体表面在柱体局部坐标系中的径向坐标.

1.2 波场分解

根据绕射波浪理论[4]:当波浪入射后，波浪会在大尺寸结构物表面产生散射波，入射波场和散射波场叠加后，形成的新的波场.波浪入射后，多柱体阵列中柱体周围的波场不同于单柱体周围的波场，多柱体阵列中柱体上的波浪作用需要考虑入射波和散射波，其中散射波还需考虑其在柱间多次绕射，柱体上总的波浪力可以写为

式中，pSc不仅考虑了入射波浪在该柱体上引起的散射波，而且考虑了阵列中其他柱产生散射波、该散射波引起的高次散射波，因此pSc可以写为

对式(3)应用分离变量法，p(ξ，η，z) 可以写为

因此柱体总的波浪力可以写为

1.3 入射波浪压力

如图1 所示，流场中任意一点Q的入射波浪压力可以表达为

式中Xi和Yi是i号局部坐标系原点在整体坐标系中的坐标，xi和yi是Q点在i号局部坐标系中的坐标.

本文采用Abramowitz 和Stegun[28]提出的符号记法，因此在椭圆柱坐标系下，入射波浪压力可以表示为[25]

2 散射波浪压力求解

2.1 散射波浪压力

波浪入射到i号柱体后，产生的第一次散射波满足式(3)～式 (5)和式(7)，其一般表达式可以写为

将式(15)和式(16)代入式(6)，求得

同一个流场中的其他柱产生的第一次散射波也会入射到i号柱体上，表示为.所有柱产生的散射波再入射到i柱后，波浪压力可以写为

根据马蒂厄函数的加法定理[29-30]，将式(18)中j号局部坐标(ξj，μj)变换到i号局部坐标系(ξi，μi)下，如下所示.马蒂厄函数的加法定理详见附录1

2.2 高次散射波浪压力

多柱体阵列中低次散射波在柱间传播，并在柱体表面产生高次散射波.其他柱产生的q-1 次散射波使i号柱体产生第q次散射波，i号柱体q次散射波浪力可以表示为

如式(18)，流场中其他柱产生的q-1 次散射波传播到i号柱体上产生的波浪压力可以写为

将式(21)转换坐标到i柱局部坐标系下为

为求式(20)中的待定系数，将式(20)和式(22)代入式(6)，得到待定系数如下

2.3 总波浪力

多柱体体系中，i号柱体上总的波浪压力可以写为

i号柱体所受的总波浪力可以表示为

3 数值算例

3.1 方法验证

首先，由于本方法在计算求解时会对式(25)中高次散射波的累加进行截断，而截断会引入一定的截断误差，所以先讨论该截断误差对计算结果的影响;其次，通过Wang 等[14]提出的波浪压力数值解验证本文提出的椭圆形截面柱体波浪压力的解析解.四柱体阵列如图2 所示，图2 中D为柱体间净距，定义其与椭圆长轴的比值Dr=1.0，4 个椭圆柱柱体具有相同的尺寸，长轴和短轴比值a/b=1.5.图3 给出了图2 布置下不同截断数q时，随波数ka变化的曲线，可以看出当q＞ 12 时本方法获得波浪力计算结果基本稳定，因此为使解析解求解高效且准确，后续计算采用q=15 为高次散射波的截断数.图4 为ka=1时，本文解析方法和数值方法得到的柱体上波浪压力，图5 为Dr=1 时，本文解析方法和数值方法得到的柱体上波浪力，假设波浪以两个入射角α入射;图6 为ka=1 时，本文解析方法和数值方法得到的流场波浪压力云图，图7 为解析方法与数值方法计算得到的C1 和C2 柱上波浪压力的相对误差;可以看出本文解和椭圆数值解相对误差较小，各图中本文解结果和数值解结果吻合良好.

图2 四柱体阵列图Fig.2 Sketch of four bodies arranged in a square form.

图3 截断误差对本文解计算结果的影响Fig.3 Impacts of truncation error on the present method

图4 不同入射角(α=0° 和 α=90°)下数值解[14]与本文解的对比Fig.4 The wave pressures on bodies Piversus θ with N=4 and ka=1 obtained present method by and the FEM[14]

图5 不同入射角(α=0° 和 α=90°)下数值解[14]与本文解的对比Fig.5 The wave forces on bodies Fi versus ka with N=4 obtained by the present method and FEM[14]

图6 入射角 α=0° 时数值解[14]与本文解云图Fig.6 The wave fields with N=4，α=0° and ka=1 obtained by the FEM[14] and present method

图7 入射角 α=0° 时数值解与本文解的相对误差Fig.7 The relative error of the wave pressure between the present method and the FEM with N=4，α=0° and ka=1

3.2 方法应用

本节将讨论高次波的影响，根据柱体排列设置了两种柱体阵列:①双柱(图8)、②四柱(图2).所示各工况中柱体截面尺寸相同，长轴和短轴之比a/b=1.5.q=2 代表的曲线是不考虑高次散射波的解析解，q=15 代表的曲线是考虑高次散射波的解析解.从图9 可以观察到，单柱情况下的波浪力和多柱情况下的波浪力是不同的，该现象随着柱数的增加而显著，且和椭圆长轴与水流方向的夹角有关.

图8 双柱体阵列图Fig.8 Sketch of arrangement of twin bodies standing side by side

图10 给出了双柱阵列中，高次散射波对柱体上波浪作用影响的结果.图中纵坐标为波浪力比值表示各柱体上考虑高次散射波的总波浪力，表示各柱体上不考虑高次散射波的总波浪力.图中给出了3 种波浪入射角度下(0°，45°和90°)，ka取值在0.2～2.0 之间的结果.图9(a)中，因为阵列沿x轴对称所以C1 和C2 上的波浪力比值相等，高次波的影响不明显;从图10 (a)和图10(b)中可以看出，当ka小于0.5 时，高次波对两柱体的影响可以忽略，当ka大于0.5 后，高次波的影响增大到不可忽略.图11 给出了ka=1 时Dr取值在0.5～5.0 之间的结果.图11(b) 和图11(c) 表明Dr大于2 后，二柱阵列的柱体受高次波的影响才会趋于减少，但即使净距为柱体长轴两倍以上，高次波的影响仍然明显;图11(c)中，C1 受到上游柱体C2 的保护，C1 受到高次波的影响比C2 小.

图9 不同入射角(α=0° 和 α=90°)下两种阵列的本文解与单柱解的对比Fig.9 The wave forces on bodies Fi versus ka with N=2 and 4 compared with that of an isolated body

图10 双柱阵列中各柱体总受力比值Fig.10 Scaling values of total wave force on bodies versus ka as twin bodies standing side by side

图11 双柱阵列中各柱体总受力比值Fig.11 Scaling values of total wave force on bodies versus Dr as twin bodies standing side by side

图12 给出了四柱阵列中，高次散射波对柱体上波浪作用影响的结果.由于对称性，波浪沿x轴入射时(α=0°)，C1 和C3 计算结果相同，C2 和C4 计算结果相同;同理，波浪沿y轴入射时(α=90°)，C1 和C2 计算结果相同，C3 和C4 计算结果相同.比较两个入射波浪方向不同的结果，可以发现因为柱体在x轴上的投影面积大于柱体在y轴上的投影面积，所以波浪沿y轴入射时高次波对柱体所受波浪力的影响更大.图13 给出了ka=1 时Dr取值在0.5～5.0 之间的结果.四柱阵列的下游柱体波浪力比值在0.5 ＜Dr＜ 1.0 间大于上游柱体的比值，在Dr＞ 2.0 之后才小于上游柱体的比值并趋于减小.

图12 四柱阵列中各柱体总受力比值Fig.12 Scaling values of total wave force on bodies versus ka as four bodies arranged in a square form

图13 四柱阵列中各柱体总受力比值Fig.13 Scaling values of total wave force on bodies versus Dr as four bodies arranged in a square form

图14 和图15 给出了柱体数量不同的情况下，C1 上波浪力比值随参数(ka，Dr)的变化.综合来看，随ka和Dr增加，四柱阵列中C1 波浪力比值变化，比双柱阵列情况下的比值变化要剧烈;随着柱体数量的增多，由于叠加效应，高次波影响也随之增加.

图14 两种阵列中C1 柱体总受力比值Fig.14 Scaling values of C1 versus ka in two arrangements

图15 两种阵列中C1 柱体总受力比值Fig.15 Scaling values of C1 versus Dr in two arrangements

图15 两种阵列中C1 柱体总受力比值(续)Fig.15 Scaling values of C1 versus Dr in two arrangements (continued)

表1 对比了本文解和参考解[14]在ka=1 时，不同柱数，不同柱间距时波浪力计算效率.因为有限元方法计算耗时和网格单元划分大小以及网格数目成正相关，所以当柱数增加或者柱间距增加后，用来模拟水域的单元也随之增多，所以本文解比参考解效率高.

表1 计算波浪力的用时 (s)Table 1 The numerical costs for calculating the total wave forces (s)

4 结论

根据绕射波理论等，基于椭圆柱坐标系，首先通过求解马蒂厄方程，得到椭圆单柱体结构波浪压力公式，再考虑多柱体体系中高次散射波问题，推导得到多柱体体系中椭圆柱体结构波浪力计算公式.本文方法与已有数值方法对比结果表明，本文解和数值解吻合的较好，而当柱数目增加或柱间距增大时，本文解的计算效率比有限元方法的高.

将本文方法应用于计算双柱阵列和四柱阵列的波浪力，分析了高次散射波对结构所受波浪作用的影响.结果表明:波数ka＜0.5 时，高次散射波影响较小，大波数的情况下，不能忽略高次波的影响;随着柱体间距离的增加，高次波的影响有减小的趋势，但仍存在波动;高次波对上游柱体波浪力的影响比下游柱体大;多柱体体系中，柱体数量增加后，柱体产生的高次散射波会叠加，高次波影响也随之增加，而结构所受的高次波作用因参数发生的波动会变剧烈.

附录1

Chatjigeorgiou 等[17]在Særmar k 等[21]的基础上化简了将马蒂厄函数的加法定理，具体可以表示为

式中，Jm(·) 为第一 类m阶贝塞尔函数，Ym(·) 为第二类m阶贝塞尔函数，为第一 类m阶汉克尔函数，为第二类m阶汉克尔函数.当n-p和s-m为奇数时，d′和d为0.