波浪对双层圆弧型贯底式开孔介质防波堤的绕射

林皋，刘俊

(1.大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连 116024;2.大连理工大学 建设工程学部 水利工程学院，辽宁 大连116024)

20世纪60年代由加拿大学者Jarlan率先提出了开孔板式防波堤［1］，该结构可以有效地降低结构波浪力和反射系数的值.随后许多研究者提出并研究了许多改进型的开孔防波堤结构，以满足实际工程的不同需要［2-3］.与此同时，为了在更广泛的频率范围内降低结构波浪力和反射系数的值，两层或多层开孔板式防波堤是较佳的选择.如意大利波尔图Torres海港中修建了三层开孔板沉箱结构［4］，最近，中国大连市化工港正在修建五层开孔板沉箱结构［3］.在理论研究方面，Sawaragi等［5］研究了波浪与双层开孔板沉箱作用条件下的衰减，结果表明，在特定的条件下，双层开孔板沉箱对波浪能量的吸收比单层开孔板要好.Kondo［6］发展了一种线性长波条件下的解析模型.李玉成等［7］应用特征函数展开法分析了双层局部开孔板沉箱对波浪反射.汪宏等［8］通过实验给出了双层开孔直立式消浪板不同开孔形式在规则波作用下的透射系数.刘勇等［9］利用特征函数展开方法研究了无穷层局部开孔板沉箱对波浪反射.有关更多的波浪与双层或者多层开孔结构的相互作用问题研究可详见文献［3］.但是，以上文献大多局限于二维平面波问题.程建生等［10］用特征函数展开解析方式研究过平面波作用下开孔介质圆弧贯底式防波堤的防浪效果的三维问题.然而，海浪实际上非常复杂，往往需要采用短峰波来描述［11］.迄今为止未见短峰波对双层圆弧型贯底式开孔介质防波堤的绕射问题的报道.

本文将采用一种新的半解析方法－比例边界有限元方法(SBFEM)建立短峰波对双层圆弧型贯底式开孔介质防波堤的绕射问题的分析模型.在此基础上分析了短峰波波浪参数和结构的几何参数等对整个结构绕射的影响.

1 双层弧型开孔介质防波堤物理模型

如图1所示，设在均匀水深h的海域中设置有内、外环半径分别为a和b的双层圆弧形开孔介质防波堤.取坐标系O-xyz，其中，O-xy平面位于静水面上，原点在圆筒柱的圆心处，O-z轴垂直向上.将2个圆弧延伸构建成2个虚拟圆柱面，整个流场被划分成1个无限子域Ω3和2有限子域Ω1、Ω2.根据线性理论，无限子域Ω3中的速度总势Φ3可表示为入射势和散射势ΦSA3之和，即Φ3=+;两个有限域的速度总势分别用Φ1和Φ2表达.为方便下面公式推导，用速度势Φ统一代表Φ1、Φ2、Φ3、以及.本问题为三维问题，应用分离变量法，可分离出时间项t和空间变量z，即速度势Φ(x，y，z，t)可表达为

式中:Z(z)代表z方向的函数，ω为角频率，t为时间项.满足边界条件的解为(k为波数):

图1 双层圆弧型开孔介质防波堤及比例边界有限元模型Fig.1 Double-layered arc-shaped porous breakwater and model of SBFEM

与Φ类似，无限子域Ω3中，速度总势φ3表示为入射势和散射势之和，即φ3=+有限子域可直接应用速度总势φ1、φ2.其中，、φ1和φ2均可由二维Helmholtz方程进行描述可统一地表示为φ(x，y)，即:

两有限域界面处的边界条件可由下式表达:

其中，对于多孔介质圆弧段标量G1等于其自身的孔隙影响系数，对非多孔介质圆弧段G值取为无穷大.

有限域Ω2和无限域Ω3界面处的边界条件可由下式表达:

类似地，G2为有限域Ω和无限域Ω3界面处的孔隙影响系数.

入射速度势可由下式表达:

总场速度势φ求出之后，其他变量如波速度，波面高度和波动压力等可通过以下公式求出:

2 比例边界有限元方程推导

通过以上阐述可以发现，问题归结为求解有限域的速度总势φ1、φ2和无限域的散射势，使之满足方程和相应的边界条件.仍用φ统一表示φ1、φ2和，并设φ在各子域边界上的法向导数为

方程(2)和(8)的等效泛函变分可表达为

图2 比例边界有限元方法计算Fig.2 The coordinate definition of SBFEM

应用SBFEM求解泛函变分时，首先需要建立笛卡尔坐标系(x，y)和比例边界有限元坐标系(ξ，s)之间的转换关系.在域内选择一点О作为SBFEM的相似中心(通过此中心，该域的全部边界可见，也就是说只有半径为b圆边界需要离散)，分别建立一个围绕比例中心的环向坐标s(s在数量上可表示为从E点开始逆时针旋转半径为b的圆弧长度)和从比例中心出发向计算区域边界辐射的辐射坐标ξ，如图2所示.环向坐标的定义从侧边面x=sb开始，逆时针旋转到侧边面s=sa，当sa和sb重合时，计算域为封闭的圆或圆环，本文中sa和sb是重合的.辐射坐标在比例中心处为ξ=0，在计算区域的半径为b边界上ξ=1.结合图1和图2，对于本文中计算情况:有限域Ω1中0=ξ0≤ξ=ξ1=a/b;有限域Ω2中a/b=ξ0≤ξ=ξ1≤1;无限域Ω2中1=ξ0≤ξ=ξ1=∞.在笛卡尔坐标系中计算域内的一点可表示为(^x，^y)，在边界处(ξ=1)上则表示为(x，y).相似中心处表示坐标为(x0，y0).图2中，阴影部分代表一SBFEM单元.两坐标系统的变换关系如下(s可以通过边界节点插值得出，本文采用三节点二次曲线插值函数N(s)，即s=N(s)s，s，为三节点弧长):

两坐标系统的求导，通过雅克比矩阵加以联系:

其中，

在比例边界坐标下，梯度算子可以表示为

其中:

根据等参变换概念，速度势函数可采用同样的插值函数进行离散:

将式(13)、(14)代入式(9)，经整理后可得:

其中:

引入以下系数矩阵E0、E1、E2、M0和Fs(ξ):

式中，Fs(ξ)为侧面边界(s=sa、s=sb)第二类边界条件不为零时在s=sa、s=sb处产生的通量.当= 0时或s=sb、s=sb存在第一类边界条件时，Fs(ξ)= 0，注意到本文中Fs(ξ)=0.

对式中含δφ(ξ)，ξ的项做分部积分，考虑到δφ (ξ)任意性，并令ζ=kbξ，最终可得如下简化方程［13］:

方程为各个子域的 SBFEM的基本方程，式(17)、(18)为其内、外边界条件.

3 比例边界有限元方程求解

3.1 对无限域Ω3的求解

方程为贝塞尔(Bessel)微分方程组，考虑到无穷远的 Sommerfeld辐射边界条件，其解可以Hri(ζ)Ti作为基函数来表示:

其中，rj为汉克尔函数的阶数，Tj为n阶向量，ci为参与系数，Hrj(ζ)为第一类汉克尔(Hankel)函数.

将式(19)代入方程(16)，并且考虑到汉克尔函数的导数性质，经整理后可得［13］:

式(21)为一个特征值问题.设 λj为矩阵的特征值，则，并且T为的特征向量.

由于汉克尔函数在无穷远处自动满足Sommerfeld辐射边界条件，所以只需考虑比例边界有限元在ξ=1处的边界条件:

3.2 对有限域Ω1、Ω2的求解及其和无限域Ω2耦合

类似地，有限域Ω1中的总场φ1可用Bessel函数作为解的基函数:

其中，

同时考虑SBFEM的边界条件可得:

其中:

同理有限域Ω2中的总场φ1可表示为:

考虑内边界条件(17)、外边界条件(18)分别可得Ω2中的SBFEM内外边界条件:

其中:

考虑到有限域Ω1和Ω2之间的边界条件可得:

其中:

其中，G1为多孔介质的孔隙影响系数对角矩阵，对角元素为SBFEM单元对应的孔隙影响系数，对于多孔介质圆弧段选取自身的孔隙影响系数，而对非多孔介质圆弧段选为无穷大.

并注意到有限域Ω2和无限域Ω3界面存在如下关系:

再考虑方程限域Ω2和无限域Ω3界面的协调边界条件(4)可解出C3:

其中，对角矩阵G2定义类似T1，¯vI3n为Ω3中边界的入射波节点速度，而且:

将式(32)代入式(28)可得C2，进一步将C2C3代入方程(31)、(29)可分别求出C和C1.将系数向量C、C1、C2和C3代入各自的域的基本解，可得到各自域的速度总势.

4 计算结果及分析

为了检验SBFEM的准确性和高效性.选取外接开孔结构物为闭合的圆形作为算例.该问题对入射波为平面波(ky=0)时有解析解.选取计算条件为: h=15 m，A=0.5 m，b=10 m，a/b=0.2，G=1.0，kb≈0.8.对应于本模型的参数为:双层圆弧变成为全圆、G1=0(相当不透水壁)和G2=1.图3给出了SBFEM和解析解［13］沿内外柱壁归一化波浪爬升值(|η|/A)的对比，结果表明SBFEM只需在外壁离散4个3节点单元就能很好的和解析解吻合.本文采用Intel Core Q8200(2.33GHz)计算机对算例进行了计算时间测试.即使在离散16个SBFEM单元的情况下，其消耗的时间不到1s，由此可以看出计算效率非常高.

图4进一步给出了以上5种不同角度入射波(取值分别为β=0、π/8、π/4、3π/8、π/2.β定义如下:β=arctan(ky/kx).其中，β=0对应为平面波; β=π/2则为驻波)作用下，掩护区中心点处(r=0)的归一化的波浪爬升(|η/a|)的变化规律.模拟计算中的参数为:两圆弧段的中心点在D，圆弧段的角度θ=120°(见图1)，h=15 m，A=0.5 m，b=2 m，a/b=0.5，G1=0.1，G2=1.0.从图4可以看出，波浪爬升呈现震荡趋势.随着β不断增大，掩护区域内中心处的波浪爬升增大.

图5进一步给出短峰波(kx=ky=sprt(2)/2)作用下，3种圆弧型开孔型防波堤结构(θ=120°):1)存在双层圆弧开孔结构;2)只存在单独内层圆弧开孔结构;3)只存在单独内层圆弧开孔结构，中心点处(r=0)的归一化的波浪爬升(|η/A|)对比结果.从图中可以明显发现，双层开孔结构比单层开孔结构的在更广泛的频率范围内具有更好的防护效果，同时可以发现，在中间频率段(大致ka=1.6～3.7)，单独内层和单独外层情况的中心点处的波浪爬高非常接近.

图3 平面波作用下的内、外柱壁波浪爬升Fig.3 Run-up of an incident plane wave around the interior and exterior cylinders.

图6分别表明了内、外圆弧段孔隙影响系数G1(G2=1.0)、G2(G1=1.0)变化对中心点处(R=0)的归一化的波浪爬升(|η/A|)的影响.除G1和G2变化外，其他计算参数和图5一致.从图6(a)可以看出，对随着孔隙影响系数的逐渐增大，中心点处的波浪爬升先是快速下降到最小值后，再逐步缓慢增大并最终趋于稳定值.而且随着β(除驻波外)的增大，波浪爬升最小值点往更小G2的靠近.从图6(b)可以看出，中心点处的波浪爬升的变化规律可以分成2种情况:1)对于较小β值(β=0，π/8)的波，随着G2的增大，中心点处的波浪爬升也是先减小再增大并趋于一渐进值;2)对于较β大值(β≥π/4)的波，随着G2的增大，中心点处的波浪爬升快速增大并趋于一渐近值.

图4 随相对波数变化下中心点处(r=0)波浪爬升Fig.4 Variation of wave run-up at r=0 versus

图5 存在不同结构层下的中心点处(r=0)波浪爬升对比Fig.5 Comparison of wave run-up at r=0 between different structures

图6 随孔隙影响系数G1、G2变化下中心点处波浪爬升Fig.6 Variation of wave run-up versus G1，G2

在以上两种分析结果基础上，图7给出了不同波对域内绕射波轮廓的影响.其中，分别选取3种波，即:平面波(kxa=1=1.0，kya=0)，短峰波(kxa=0，kya=1.0)，并且选择了3种不同的孔隙影响系数G1=G2分别为0.1、1.0和10.0.其他计算参数和图4一致.从图中可以看出，所有的绕射波轮廓图都关于x轴对称;随着G1和G2的减小，外圆弧段和内筒间的波浪高度呈现明显减小趋势.垂直入射平面波能体现出明显的反射特性，即掩护效果更好，但当孔隙影响系数G1和G2很大情况下，也表现得不是很明显;而驻波和短峰波的绕射波相对平面波的更为复杂.

图8表明了内外圆弧半径比a/b变化对中心点处(r=0)的归一化的波浪爬升(|η/A|)的影响.其他计算参数和图4一致.从图中可以看出，随着归一化相对波数kb/2的增大，中心点处的波浪爬升先是快速下降，降到到一定值后又呈现震荡变化.在低频区，半径比a/b对中心点处的波浪爬升影响非常小，而在中高频率段，半径比a/b只影响震荡的波峰及波谷出现的位置.

图9表明了3种组合圆弧张开角度θ变化对中心点处(r=0)归一化的波浪爬升(|η|/A)的影响.除圆弧张开角度θ外，其他计算参数和图4一致(两圆弧段的中心点在D).3种组合分别为:1)内圆弧段的角度变化而外圆弧段角度固定为120°;2)外圆弧段的角度变化而内圆弧段角度固定为120°;3)内、外圆弧段的角度同时变化.从图中可以看出，总体来说，不管那种组合情况，随着张开角度的增大，中心点处的波浪爬升震荡的振幅增大;组合情况2)的震荡的频率更大;而且可以进一步证实，组合情况3)在更大的频率范围内掩护效果更好.

图10表明了3种组合圆弧位置变化γ对对中心点处(r=0)的归一化的波浪爬升(|η|/A)的影响.圆弧中心点从A点逆时针转到D点(总共180°).其他计算参数和图6一致.3种组合分别为:1)内圆弧段的位置变化而外圆弧段位置不变;2)外圆弧段的位置变化而内圆弧段位置不变;3)内、外圆弧段的位置同时变化.从图中可以看出，不管哪种组合情况，随着r的增大，中心点处的波浪爬升震荡的振幅减小，即r值越大掩护效果更佳;组合情况2)的震荡的频率更大，同时组合情况3)的振幅更大.

图7 不同短峰波在不同孔隙影响系数绕射波的波幅分布Fig.7 Diffracted wave amplitude distributions with different incident wave and different porous-effect parameters

图8 随内、外半径比变化下中心点处波浪爬升Fig.8 Variation of wave run-up versus a/b

图9 随圆弧段张开角度变化下中心点处波浪爬升Fig.9 Variation of wave run-up at r=0 versus opening angle

图10 随圆弧段位置变化下中心点处波浪爬升Fig.10 Variation of wave run-up at r=0 versus location of arcs

5 结论

开展了应用比例边界有限元法探讨短峰波对双层圆弧型贯底式开孔介质防波堤的绕射问题.将开孔圆弧段延伸形成虚拟圆弧段，以不同的孔隙影响系数反映实际结构段与虚拟圆弧段的不同特性，从而构建出新的匹配边界条件，能使该问题落入SBFEM求解范围之内.基于比例边界坐标变换和变分原理推导出SBFEM表示的本问题的基本方程.在求解过程中，针对有限域和无限域分别引入了系列贝塞尔函数和汉克尔函数作为其解的基函数，并且结合两个有限域界面处以及有限域和无限域界面处的边界条件求解出对应的各系数向量.数值算例与解析解结果对比表明该方法具有很高的精度和计算效率.最后，详细探讨出诸如相对波数、半径比、圆弧段的位置、张角以及孔隙影响系数G对整个结构绕射的影响.结果表明:1)平面波的防浪效果最好，而驻波的效果最差;2)在更广泛的频率范围内，双层开孔结构对波浪的防浪效果比单层开孔结构要好;3)在内层的孔隙影响系数较小的情况下，内层的孔隙影响系数增大将能有效地降低中心点处的波浪爬升;4)而外层的孔隙影响系数增大将会使中心点处的波浪爬升提高;5)在较低的频率范围内，内、外柱半径比对波浪的防浪效果影响较小;6)内、外柱的张角变化以及它们的位置变化对中心点处波浪爬升影响非常大.本研究可为工程上设计双层圆弧型贯底式开孔介质防波堤提供理论参考.

［1］JARLAN G E.A perforated vertical wall breakwater［J］.The Dock and Harbor Authority，1961，41(486):394-398.

［2］CHWANG A T，CHAN A T.Interaction between porous media and wave motion［J］.Annual Review of Fluid Mechanics，1998，30:53-84.

［3］HUANG Z H，LI Y C，LIU Y.Hydraulic performance and wave loadings of perforated/slotted coastal structures:a review［J］.Ocean Engineering，2011，38(10):1031-1053.

［4］FRANCE L.Vertical breakwaters:the Italian experience［J］.Coast Engineering，1994，22:31-55

［5］SAWARAGI T，IWATA K.Wave attenuation of a vertical breakwater with two air chambers［J］.Coastal Engineering in Japan，1978，21:63-74.

［6］KONDO H.Analysis of breakwaters having two porous walls［C］//Proceedings of Coastal Structures’79.Virginia，USA，1979.

［7］李玉成，刘洪杰，滕斌.双层局部开孔板沉箱对波浪反射的理论研究［J］.海洋工程，2005，23(1):18-27.

LI Yucheng，LIU Hongjie，TENG Bin.Theoretical analysis of the wave reflection by caissons with double-layered perforated wall［J］.The Ocean Engineering，2005，23(1):18-27.

［8］汪宏，沈丽玉，王勇.双层开孔直立式板结构的消波性能试验［J］.水运工程，2011，450(2):21-25.

WANG Hong，SHEN Liyu，WANG Yong.Test on wavedissipating performance of double-layer-hole vertical plate structure［J］.Port＆Waterway Engineering，2011，450 (2):21-25.

［9］LIU Y，LI Y C，TENG B.Interaction between obliquely incident waves and an infinite array of multi-chamber perforated caissons［J］.Journal of Engineering Mathematics，DOI 10.1007/s10665-011-9484-2.

［10］程建生，缪国平，王景全，等.圆弧型贯底式多孔介质防波堤防浪效果的解析研究［J］.哈尔滨工程大学学报，2008，29(1):5-10.

CHENG Jiansheng，MIAO Guoping，WANG Jingquan，et al Wave-prevention efficiency of an arc-shaped bottom-mounted porous breakwater［J］.Journal of Harbin Engineering University，2008，29(1):5-10.

［11］TSAI C P，JENG D S，HSU J R C.Computations of the almost highest short-crested waves in deep water［J］.Applied Ocean Research，1994，16(6):317-326.

［12］WOLF J P，SONG C M.Consistent infinitesimal finite element cell method:three-dimensional vector wave equation［J］.International Journal for Numerical Methods in Engineering，1996，39(13):2189-2208.

［13］TAO L B，SONG H，CHAKRABARTI S.Scaled boundary FEM model for interaction of short-crested waves with a concentric porous cylindrical structure［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，2009，135 (5):200-212.