短周期规则波浪在透水鱼礁型潜堤群上的传播特性

摘要

为探究系列以周期性分布的单排透水鱼礁型潜堤群对波浪传播变形以及波生湍流场的响应特性，通过水槽试验和数值模拟，对短周期规则波作用下不同间距布设的鱼礁型潜堤群波浪传播特性展开了研究。试验测定了短周期入射波况下无间距鱼礁型潜堤群的波面高度及其变化过程，波浪在发生浅水变形后波能耗散波高显著减小，但由于堤顶淹没水深较大波浪仍保持准对称性的正弦波形态向前传播。基于光滑粒子流体动力学（smoothed particle hydrodynamic，SPH）方法建立了验证后的二维无反射数值波浪水槽，模拟了间距比为0.5~3.0的鱼礁型潜堤群的波浪传播过程。数值结果表明：鱼礁型潜堤群的间距对波浪反射影响非常大，各工况下鱼礁型潜堤群的反射系数均大于0.75；透射系数总体上随着间距比的增大而减小。特别地，当间距比为2.5时，堤顶波高骤升，反射系数达到最大，透射系数几乎最小，这是由于波浪与鱼礁型潜堤群发生了同频共振的布拉格反射效应，此时迎浪堤上方的上涌流范围和流速也达到最大。

关键词

鱼礁型潜堤; 光滑粒子流体动力学; 布拉格反射; 波生流

近年来，生态学在海岸工程中日益受到重视，效法自然的工程设计理念成为主流。鱼礁型潜堤是一类新兴的生态型海岸防护构筑物，将人工鱼礁的设计融合到潜堤的结构中^［1］，兼具生态修复和消减波浪的效果。2010年，我国北戴河首次采用鱼礁型潜堤代替传统堆石离岸堤用于海滩恢复治理，工程后近岸海域生态环境显著改善^［2-3］。鱼礁型潜堤作为外海波浪向岸传播入湾的第一道屏障，抵御波浪的能力对当地水动力和水生态环境的塑造至关重要。

波浪在鱼礁型潜堤上的传播是一个复杂又存在一定内在规律的过程，包括越浪、波浪反射以及波浪透射等物理现象，自20世纪70年代以来一直是国内外学者研究的重点。Armono和Swamidas^［4］、Garcia等^［5］、Chen等^［6］模拟研究了单个和多个透水鱼礁型潜堤周围波流变化，均发现潜堤后方的水位有所抬升。张志涛等人^［7］通过水槽试验证明以不同间距布设的鱼礁型潜堤可以降低波浪能量并使其提前破碎。李欣雨^［8］和潘冬冬等^［9］分别构建三维和二维数值波浪模型分析了规则波作用下潜堤的流场变化，得到波高和潜堤间距对周围流场的形成均有较大的影响。Saprykina等^［10］和Metallinos等^［11］利用室内试验和数值模拟研究了放置在海底斜坡上的两潜堤在不同间距时对波浪传播的影响，结果表明双潜堤的布置能有效降低波高和平均波周期，而平均波周期的减小程度取决于潜堤间距。目前许多关于鱼礁型潜堤布设间距对波浪传播变化影响的研究大多集中在两单体之间，对3个以上组成的大规模鱼礁型潜堤群的波浪响应特征讨论较少。在最新的相关发表成果中，舒安平等^［12］分析了渤海湾不同布设间距下大规模鱼礁型潜堤群的流场分布特征，但其研究并未考虑海湾中主要动力因素——波浪。因此，在实际工程中，对于波浪在大规模鱼礁型潜堤群上传播特性的分析迫切地需要更多研究成果以提供科学理论支持。

目前，对于波浪在鱼礁型潜堤上传播过程的数值模拟以求解欧拉网格划分的N-S模型为主，在计算过程中需要不断使用网格重构技术以避免网格畸变导致的计算中断。随着现代计算机能力的提高和波浪模型算法的发展，无网格的光滑粒子动力学（smoothed particle hydrodynamics，SPH）方法在计算流体动力学方面的应用也得到了快速发展^［13-14］。SPH 模型中粒子具有自适应性，其在处理大变形自由表面流问题和与复杂结构物相互作用的物面边界处理时，显著优于欧拉网格法^［15-17］。需要注意的是，相比网格类方法，SPH方法在计算过程中计算量大、耗时长，近年来利用GPU（graphical processing unit）并行技术加速SPH模型的算法开发已大大改善了这种情况^［18-19］，但对于工程尺度下的数值模型计算，SPH方法仍然面临一定挑战。

本文以北戴河生态修复工程中投放的透水鱼礁型潜堤为原型，通过合理比尺（1：10）缩放制作试验模型，对由多个单体潜堤模型组合成的大尺度鱼礁型潜堤群开展物理模型试验，讨论分析了无间距鱼礁型潜堤群上的波高空间分布和结构消波特性；基于DualSPHysic开源平台，采用SPH数值方法，建立与物理模型试验相同的数值波浪模型，通过改变鱼礁型潜堤群的布设间距，模拟分析不同间距比

（d/l_r=0.5、1.0、1.5、2.0、2.5、3.0）下的鱼礁型潜堤群对规则波浪的响应特征和波生紊流场的反馈机制，为进一步了解和优化以周期性间距分布的系列鱼礁型潜堤群在海岸生态修复工程中的应用提供科学依据。

1 SPH数学模型

SPH是以拉格朗日视角出发的无网格数值计算方法，将整个连续的计算域离散为N个具有质量、密度和位置等物理属性的粒子。在计算过程中，根据周围粒子的物理性质，通过离散形式的Navier-Stokes方程（简称N-S方程）对每个粒子的位置进行局部积分，随着相邻粒子间的相互作用，每一时刻的流体粒子位置都在不断更新。

1.1　水动力控制方程

本文主要研究重力水波的运动过程，在SPH方法中水体被视作微可压缩流体，基于N-S方程，SPH中流体的质量守恒和动量守恒方程可以表达为

$$\frac{\mathrm{d}{\rho }_{\mathrm{a}}}{\mathrm{d}t}={\rho }_{\mathrm{a}}\sum _{\mathrm{b}}{m}_{\mathrm{b}}{v}_{\mathrm{a}\mathrm{b}}{\nabla }_{\mathrm{a}}{W}_{\mathrm{a}\mathrm{b}}$$

（1）

$$\begin{array}{l}\frac{\mathrm{d}{v}_{\mathrm{a}}}{\mathrm{d}t}=\\ -\sum _{\mathrm{b}}{m}_{\mathrm{b}}(\frac{P_{\mathrm{b}}}{{\rho }_{\mathrm{b}}^{\mathrm{2}}}+\frac{P_{\mathrm{a}}}{{\rho }_{\mathrm{a}}^{\mathrm{2}}}+{\mathrm{\Pi }}_{\mathrm{a}\mathrm{b}}){\nabla }_{\mathrm{a}}{W}_{\mathrm{a}\mathrm{b}}+g\end{array}$$

（2）

式中：下标a代表任意流体粒子，下标b代表a粒子的邻域粒子，它包括流体粒子与边界粒子；ρ_a、m_a、v_a与P_a分别是任意流体粒子a的密度、质量、速度与压力；v_ab表示粒子a与b之间的相对速度；∇_a表示相对粒子a坐标的梯度；W_ab为粒子a对粒子b产生影响的光滑核函数，W_ab=W（r_a-r_b， h），r为粒子的空间位置，h为光滑长度。

当考虑流体黏性耗散时，式（2）中的人工黏度项表达^［20］为

$${\prod }_{\mathrm{a}\mathrm{b}}=\left\{\begin{array}{l}\frac{-\alpha \overline{c_{\mathrm{a}\mathrm{b}}}{\mu }_{\mathrm{a}\mathrm{b}}}{\overline{{\rho }_{\mathrm{a}\mathrm{b}}}},v_{\mathrm{a}\mathrm{b}}{r}_{\mathrm{a}\mathrm{b}}<\mathrm{0}\\ \mathrm{0},v_{\mathrm{a}\mathrm{b}}{r}_{\mathrm{a}\mathrm{b}}>\mathrm{0}\end{array}\right.$$

（3）

式中：$\overline{c_{\mathrm{a}\mathrm{b}}}$表示粒子a与b声速的算术平均值；α为调整数值耗散的经验常数，对于研究波浪在沿海结构物上传播和波浪荷载的问题，α取0.01；为避免粒子相互靠近时导致数值发散引入的计算因子，${\mu }_{\mathrm{a}\mathrm{b}}=h{v}_{\mathrm{a}\mathrm{b}}\cdot \frac{r_{\mathrm{a}\mathrm{b}}}{{r_{\mathrm{a}\mathrm{b}}}^{\mathrm{2}}+{\phi }^{\mathrm{2}}}$，$\phi =\mathrm{0.1}h$。

流体压力采用Monaghan^［21］提出的计算自由表面流动的状态方程如下：

$$P=B\left[{\left(\frac{\rho }{{\rho }_{\mathrm{0}}}\right)}^{\mathrm{7}}-\mathrm{1}\right]$$

（4）

式中：P为流体初始密度，本文取水体密度1 000kg·m^-3；B为限制密度场的最大改变量，参考声速 的值要求至少比最大流体速度大10倍，以控制流体密度变化率在1%以内^［21］。

基于上述SPH形式的质量和动量方程，为得到更稳定的压力场和速度场，经泰勒展开修正后的质量和动量守恒方程为

$$\frac{\mathrm{d}{\rho }_{\mathrm{a}}}{\mathrm{d}t}=\frac{{\rho }_{\mathrm{a}}{\displaystyle \sum _{\mathrm{b}=\mathrm{1}}^N}\frac{{m^{\mathrm{*}}}_{\mathrm{a}\mathrm{b}}}{{{\rho }^{\mathrm{*}}}_{\mathrm{a}\mathrm{b}}}\left({u^R}_{\mathrm{a}}-{u^R}_{\mathrm{a}\mathrm{b}}\right){\nabla }_{\mathrm{a}}{W}_{\mathrm{a}\mathrm{b}}}{{\displaystyle \sum _{\mathrm{b}=\mathrm{1}}^N}\left({r^{\mathrm{*}}}_{\mathrm{a}\mathrm{b}}-r_{\mathrm{a}\mathrm{b}}\right)\otimes {\nabla }_{\mathrm{a}}{W}_{\mathrm{a}\mathrm{b}}\frac{{m^{\mathrm{*}}}_{\mathrm{a}\mathrm{b}}}{{{\rho }^{\mathrm{*}}}_{\mathrm{a}\mathrm{b}}}}$$

（5）

$$\frac{\mathrm{d}{u}_a}{\mathrm{d}t}=\frac{-{\displaystyle \sum _{\mathrm{b}=\mathrm{1}}^N}{m^{\mathrm{*}}}_{\mathrm{a}\mathrm{b}}\frac{{P^{\mathrm{*}}}_{\mathrm{a}\mathrm{b}}}{{{\rho }^{\mathrm{*}}}_{\mathrm{a}\mathrm{b}}{\rho }_{\mathrm{a}}}\frac{\partial W_{\mathrm{a}\mathrm{b}}}{\partial x_{\mathrm{a}}}}{{\displaystyle \sum _{b=\mathrm{1}}^N}({r^{\mathrm{*}}}_{\mathrm{a}\mathrm{b}}-r_{\mathrm{a}\mathrm{b}})\otimes {\nabla }_{\mathrm{a}}{W}_{\mathrm{a}\mathrm{b}}\frac{{m^{\mathrm{*}}}_{\mathrm{a}\mathrm{b}}}{{{\rho }^{\mathrm{*}}}_{\mathrm{a}\mathrm{b}}}}$$

（6）

$${r^{\mathrm{*}}}_{\mathrm{a}\mathrm{b}}-r_{\mathrm{a}\mathrm{b}}=\left(r_{\mathrm{b}}-r_{\mathrm{a}}\right)\frac{m_{\mathrm{a}}/{\rho }_{\mathrm{a}}}{m_{\mathrm{a}}/{\rho }_{\mathrm{a}}+m_{\mathrm{b}}/{\rho }_{\mathrm{b}}}+{u^R}_{\mathrm{a}\mathrm{b}}\mathrm{\Delta }t$$

（7）

式中：N为计算域内的总粒子数； 为张量积符号；∆t为时间步长；u_a、x_a为粒子a在x方向的速度和位置；u^R_a为粒子a的径向速度；u^R_ab为粒子a和粒子b接触点处的法向速度；P^*_ab为接触点处的压力；m^*_ab为粒子a和粒子b的算术平均值；ρ^*_ab为粒子a和粒子b密度的体积平均；${\nabla }_{\mathrm{a}}{W}_{\mathrm{a}\mathrm{b}}$为核函数对粒子a的求导（核函数梯度）。

1.2　边界条件

SPH方法是典型的拉格朗日方法，自由表面的边界条件自动满足。模型中的固壁边界条件采用动力边界方法^［22］处理，计算域内鱼礁型潜堤的壁面以及水槽的边壁和底部离散为两层固壁边界粒子。边界粒子和流体粒子一样具有质量、密度和压力等物理属性。在计算过程中，固壁边界粒子的密度和压力根据连续性方程的状态方程随时间不断更新，但其位置和速度始终保持不变，且不受作用力大小及方向的限制^［23］。动力边界法无需进行显式的边界处理，而是将边界处理耦合在控制方程的求解过程中，因此代码实现简单的同时，计算量相对较小，数值稳定性强，适用于具有复杂边界的波浪-结构物相互作用的问题。方程采用Symplectic数值积分方法求解，具有二阶精度^［24］。

数值模型采用推板造波形式，造波板由固边界粒子组成，推板造波机通过将一组控制信号输入造波机，使推波板按照一定的轨迹前后移动，从而带动水体产生波浪。其中，所生成波浪的波高和周期取决于造波板的位置时间序列，其反映了推波板同时间的位移变化关系。根据二阶造波理论，对于不可压缩、无旋的理想流体，当自由表面流体压力为常数时，波面方程和造波板位置时序列之间的传递函数为^［25］

$$Q=\frac{H}{S_{\mathrm{0}}}=\frac{\mathrm{4}\mathrm{s}\mathrm{i}\mathrm{n}{\mathrm{h}}^{\mathrm{2}}(kd)}{\mathrm{2}kd+\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{2}kd)}$$

（8）

式中：H为目标波高；S₀为推波板的运动振幅；Q代表水动力传递函数；k为波数；d为静水深。

因此，t时刻造波板位移为

$$e_{\mathrm{0}}(t)=\frac{S_{\mathrm{0}}}{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}(\omega t+\delta )+[(\frac{H^{\mathrm{2}}}{\mathrm{32}d})\cdot (\frac{\mathrm{3}\mathrm{c}\mathrm{o}\mathrm{s}h(kd)}{\mathrm{s}\mathrm{i}\mathrm{n}{h}^{\mathrm{3}}(kd)})\frac{\mathrm{2}}{m_{\mathrm{1}}}]\cdot \mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\omega t+\mathrm{2}\delta )$$

（9）

式中：ω=2π/T，T为波周期；δ为初始相位；m₁=H/S₀.

为消除波浪在造波板处二次反射的影响，缩短数值水槽长度，减少波浪传播过程中产生的波能耗散，本模型采用Madsen^［26］提出的主动吸收式数值造波技术来生成行进波。基于t时刻的造波板位置e（t），利用造波板的速度信号对其位置进行更新修正，得到修正后造波板的位置为

$$e_{\mathrm{0}}(t+\mathrm{d}t)=e_{\mathrm{0}}(t)+(U_{\mathrm{C}}(t+\mathrm{d}t)+U_{\mathrm{C}}(t))\frac{\mathrm{d}t}{\mathrm{2}}$$

（10）

式中：U_C为修正后的造波板速度，U${}_{\mathrm{C}}$（t+dt）=U_I（t）-U_R（t），U_I为造波板理论速度，U_R为吸收反射波所需的造波板修正速度；U_I（t）=ωS₀/2cos（ωt+δ），U_R（t）=（ƞ_SPH（t）-ƞ_I（t））·（g/d）^1/2；η_SPH为造波板前的实测水面高程，η_I为目标入射波的理论波面高程。

同时，在水槽末端设置一个与物模试验尺寸相同的黏性阻尼消波区，以在每个时间步结束时衰减消波区域内流体粒子的速度和压力梯度。消波区流体粒子的速度和压力表达式为

$$v_{\mathrm{a}}=v_{\mathrm{0}}-\mathrm{0.5}\frac{x-x_{\mathrm{0}}}{l}{\displaystyle \sum _{b=\mathrm{1}}^N}[\frac{m_{\mathrm{b}}}{{\rho }_{\mathrm{a}\mathrm{b}}}(v_b-{\nu }_{\mathrm{a}}){\displaystyle \sum _{b=\mathrm{1}}^N}{W}_{\mathrm{a}\mathrm{b}}\frac{m_{\mathrm{b}}}{{\rho }_{\mathrm{b}}}$$

（11）

$${P^{\mathrm{*}}}_{\mathrm{a}\mathrm{b}}=\frac{P_{\mathrm{a}}{\rho }_{\mathrm{b}}{c}_{\mathrm{b}}+P_{\mathrm{b}}{\rho }_{\mathrm{a}}{c}_{\mathrm{a}}+(\beta \frac{x-x_{\mathrm{0}}}{l}){\rho }_{\mathrm{a}}{c}_{\mathrm{a}}{\rho }_{\mathrm{b}}{c}_{\mathrm{b}}({u^R}_{\mathrm{b}}-{u^R}_{\mathrm{a}})}{{\rho }_{\mathrm{a}}{c}_{\mathrm{a}}+{\rho }_{\mathrm{b}}{c}_{\mathrm{b}}}$$

（12）

式中：x₀为消波区起始横坐标；x为消波区内流体粒子横坐标；l为消波区域长度。

2 试验设计

物模试验在同济大学水利与港口工程实验室的波浪水槽中进行，该水槽长50m、宽0.8m、高1.2m，最低工作水深0.2m，最高工作水深0.9m。水槽左侧安装有液压伺服电机驱动式造波机，造波系统由大连理工大学海洋工程研究所仪器设备研发组设计建造，采用基于二次谐波的主动吸收系统，造波机后侧设有消能网，以防止溅水，水槽另一端设置有消浪网，用于吸收反射波能，减弱波浪反射的影响。

试验以北戴河生态修复工程中投放的透水鱼礁型潜堤为原型进行试验设计。参考《波浪模型试验规程（JTJ/T234-2001）》相关说明，按照重力相似准则进行试验参数设计，即Fr_p = Fr_m。综合考虑目标投放海域水深、现场波浪条件、实验室设备条件和仪器测量精度等因素，确定模型几何比尺为λ_l = l_P/l_m = 10。根据前人研究成果，有机玻璃制备简易且透光率高，易于观察水力现象，适合制作室内模型^［12， 27-28］。因此，采用有机玻璃制作鱼礁型潜堤模型（单体结构尺寸为0.3m×0.2m×0.3m，表面开孔率为5%），并利用玻璃胶将多个单体潜堤模型黏合组成整体长1.8m、高0.3m的鱼礁型潜堤群，再将其用玻璃胶固定放置在距离造波机32.5m处，模型示意图及结构几何参数详见图1。

图1  波浪水槽与鱼礁型潜堤布置示意图

Fig.1  Schematic diagram of wave flume and reef-type breakwaters

波浪测量采用南京水利科学研究院研发生产的CBY-II型波高测量控制系统，试验前，先对浪高仪进行率定，保证浪高仪量程。浪高仪精度为±2 mm，采样频率为50Hz。试验共布置10个浪高测点，各测点具体位置及相互间的间距如图1所示。基于Goda两点法^［29］，在潜堤前设置2个浪高仪W_a和W_b，W_b距离潜堤迎浪面2.5m，以排除潜堤附近紊动的影响。基于波浪理论，由于潜堤位于平底上，透射波高可定义为堤后约1倍波长处的波高，故W_c设置于潜堤后5m处用于测量透射波。为测量鱼礁型潜堤周围波幅的沿程变化，浪高仪W1~W3位于潜堤前的向海侧，W4~W5设于堤顶，用于采集波浪的浅化和破碎，W6~W7位于潜堤后方的向岸侧，用于监测波浪越过潜堤后的传播变形过程。由于本研究的目的是考察以不同间距布置的鱼礁型潜堤群上波浪传播变化，因此仅选择一种代表性的水力条件（即堤顶始终淹没于水中）。实际工程中，鱼礁型潜堤投放于水深约4.7m的海域中^［30］，按几何比尺1：10将试验水深h设为0.47m，对应潜堤淹没水深h_r为0.17m。考虑水下透空结构在波浪周期较小时消波效果明显^［31］，波浪传播过程中的非线性变化剧烈，为测试鱼礁型潜堤群的消波性能，入射波周期设为T=1.2s，初始波高H₀=0.04m。水深波长比约为0.235，属于有限水深范畴。试验重复3次，取均值作为最终结果。

为研究大尺度鱼礁型潜堤群的波浪响应机制与潜堤布设间距之间的关系，引入间距比d/l_r这一参数（d为相邻两鱼礁型潜堤之间的距离，l_r为单体潜堤长度，量纲为1）。基于水槽试验结果，建立验证后的数值波浪水槽模型，如图2所示，通过改变相邻两潜堤之间的距离，模拟分析间距比d/l_r=0.5、1.0、1.5、2.0、2.5和3.0这6种工况下鱼礁型潜堤群的波浪传播变形特性，各工况下相邻潜堤前后均以等间距布置。随着潜堤布设间距的延长，测点W5~W7和W_c的位置相应向后移动，各测点相对潜堤的距离保持不变。数值实验中，模型设置的波浪周期、初始波高、水深等参数与水槽试验一致。

图2  数值波浪水槽模型及固壁边界粒子分布示意图

Fig.2  Schematic diagram of numerical wave flume model and boundary particles distribution

3 波浪在鱼礁型潜堤上的传播特性分析

3.1　波面沿程变化

在SPH模型中，初始粒子间距（d_x）表征模型分辨率，也决定了流体和固体边界粒子的总数。为准确捕捉波面运动轨迹，首先测试不同d_x对数值模型计算精度的敏感性。根据Roselli等人^［31］的研究成果，兼顾模型计算精度和时间成本，初始粒子间距应大于1/10的初始入射波高（H₀），即d_x/H₀ >1/10；更小的d_x不会再显著提高模型计算精度，甚至会使之有所下降，而计算时间会大幅增加^［32］。为便于分析粒子间距影响，引入相对粒子间距ζ这一无量纲参数，ζ=d_x/H₀（d_x为初始粒子间距，H₀为初始入射波高）。因多孔透水鱼礁型潜堤结构较为复杂，孔口、内室及固壁交界处水力条件复杂，这些区域的水力特性需要相对小的粒子间距才能更好呈现。参考前人对SPH模型初始粒子间距选取和率定情况^［15， 17， 23， 33］，考虑计算能力条件和数值波浪水槽及鱼礁型潜堤模型的结构尺寸，分别选择初始粒子间距为0.025m、0.015m和0.006m等三种情况模拟计算波浪SPH数值模型进行初始粒子间距影响分析，相对粒子间距分别为ζ=0.625、0.375和0.150。

图3为不同测点处波面历时曲线的计算值与试验值的对比图，其中计算值包含了三种粒子间距的模拟结果。从图中波面基本稳定后的历时曲线中，可以看到规则波沿鱼礁型潜堤的传播变化规律。波浪在到达鱼礁型潜堤前（W1~W3），不同初始粒子间距下模型的计算结果与试验结果均吻合较好。波浪传播至浪高仪W4~W5处水深改变，由于堤顶淹没水深较大且波浪传播较快，波浪在礁面上并没有发生明显的破碎现象，波面总体上还是呈现为正弦波的形态，波形保持对称性，但ζ=0.625下的模型计算值与试验值开始出现较大偏差。从W3和W6的波面历时曲线可以明显的观察到由鱼礁型潜堤引起的波浪反射和绕射现象，在测点W3处入射波和从潜堤迎浪面反射回来的波叠加，导致该点波高增大。对比模拟结果发现，在测点W6处由绕射形成的波浪在ζ=0.150的数值解与试验结果吻合最好，但ζ=0.625模型的计算值与实测值出现更大偏差。波浪越过潜堤后的测点W6和W7波高显著减小，说明部分波能在越堤过程中有所损耗。在测点W7处，ζ=0.375和0.150模型均与试验结果吻合较好，而ζ=0.625模型计算值与实测值之间仍然存在较大偏差。综上所述，根据波面历时变化曲线，本文建立的SPH波浪模型（ζ=0.375和0.150）能够较为准确地模拟潜堤上波高的空间分布特征，说明模型能够重现鱼礁型潜堤上波浪浅水变形的过程，同时亦能较为准确地捕捉潜堤迎浪面和背浪面的波面曲线。

图3  不同测点处波面历时曲线的SPH结果与试验结果对比

Fig.3  Comparisons of the free surface elevations between the experimental data and the SPH results at different wave gauges

为了定量评价数值和试验结果之间的误差，分别计算不同粒子间距下的数值结果与实测值间的均方根误差（RMSE），计算结果如表1所示。

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt[]{\frac{\mathrm{1}}{n}{\sum }_i^n{\left[{\eta }_{\mathrm{N}\mathrm{u}\mathrm{m}}(i)-{\eta }_{\mathrm{E}\mathrm{x}\mathrm{p}}(i)\right]}^{\mathrm{2}}}$$

（12）

式中，n为样本个数，η_Num和η_Exp分别为数值计算和试验测量得到的波高。

表1  不同粒子间距下模型计算值与实测值间误差

Tab.1  The RMSE between SPH results and experimental data for different relative spacing distance

相对粒子间距ζ	0.625	0.375	0.150
初始粒子间距dx / m	0.025	0.015	0.006
RMSE	0.013	0.009	0.005

从表1中可以看出，随着初始粒子间距的减小，模型的粒子总数和计算精度也相应提高。当ζ=0.375时，SPH模型的计算结果已经收敛；而当ζ=0.625时，由于粒子间距较大，初始粒子数量较少，导致模型预测的自由液面与试验结果偏差较大，故此粒子间距（ζ=0.625）设置不予采用。结合图3的波面变化结果，针对本文研究的鱼礁型潜堤-波浪SPH数值模型，初始相对粒子间距在0.375≥ζ≥ 0.150范围内可以准较确地模拟波浪传播变形过程。考虑到鱼礁型潜堤的结构复杂，表面有尺度较小的透水圆孔，为精确模拟结构的透水性对波浪传播过程的影响，选择ζ=0.150为模型的全局初始相对粒子间距，各工况模拟时间均为100s。

波浪的传播变形过程直接影响着滩面上的泥沙运动以及海滩剖面的塑造。由于鱼礁型潜堤的存在，部分波能在潜堤迎浪面前反射回外海，反射波与入射波叠加在迎浪面前会形成部分立波。为分析波浪在鱼礁型潜堤上的传播变形过程，以布设间距比d/l_r=0.5的鱼礁型潜堤群为例，在波浪传播稳定后，提取了一个周期内波浪的传播和变形过程，如图4所示。t=0.0T时，波浪传播至潜堤正上方，由于水深急剧变浅和潜堤的反射作用，堤前波高迅速增大，波浪明显发生浅水变形，堤顶处水体呈上涌趋势，堤前水体以较高流速通过迎浪面圆孔进入结构内部。t=0.25T时，波浪越过堤缘并在堤顶上发生变形，波高明显减小，波峰传播至堤间，波浪仍以正弦波的形态沿潜堤向前传播，此时堤顶上方形成的上涌流流速达到最大，部分上涌流通过堤顶的透水孔进入结构内部，但堤内及附近水质点仍保持较低流速。t=0.5T时，波峰变得平缓，波峰顶端流体质点运动速度降低，潜堤上方水体的运动方向与波浪传播方向保持一致，同时潜堤背浪面后出现回流现象，部分结构内部的水体通过迎浪面透水孔流出潜堤外。t=0.75T时，波浪在连续越过多个鱼礁型潜堤后波能显著损耗，波峰变得更为平缓，大部分堤顶上水体以水流的形式向前流动，潜堤内的回流紊动加剧，部分结构内的水体再次透过堤顶圆孔流出。t=T时，波浪传播进入下一个周期循环。总体来看，波浪经过鱼礁型潜堤传播的过程中，潜堤内部出现连续的回流紊动现象，使内部水体的速度值一直较小，仅为堤前速度的15%左右，相邻潜堤间的间距使波浪在传播过程中波峰水质点速度呈现增大—减小—小幅增加—再减小的变化趋势，而非单调减小，这为增强工程海域的水体交换提供了更多可能性。

图4  一个波周期内波浪经鱼礁型潜堤（d/l_r=0.5）的传播和变形过程

Fig.4  Wave propagations and deformation over reef-type breakwaters (d/l_r=0.5) in one wave period

3.2　鱼礁型潜堤对波浪反射和透射的影响

从能量角度来讲，波浪通过鱼礁型潜堤时一部分波能被反射，一部分能量在潜堤附近损耗，还有一部分能量透过潜堤。波浪反射系数K_r和透射系数K_t是评估潜堤消波性能的常用指标。在本研究给定的入射波浪条件下，波浪越过鱼礁型潜堤时未发生明显破碎现象，因此采用波高比值定义波浪的反射和透射系数。

Kr = Hr/Hi

（13）

Kt = Ht/Hi

（14）

式中：H_r为堤前反射波高；H_i为堤前入射波高；H_t为堤后透射波高。

基于能量守恒原理，如果计算反射系数的断面与计算透射系数的断面之间没有任何能量损失，那么反射系数与透射系数的平方和为1。但无论现实生活还是模型试验，这两个断面之间均存在能量损耗，由此波能耗散系数K_l的概念被提出。

Kr2+Kt2+Kl2=1

（15）

表2为不同工况下利用SPH方法模拟得到的鱼礁型潜堤反射、透射和耗散系数的计算结果。对于以一定形式布置于海底的系列鱼礁型潜堤群而言，结构体之间的间距对波浪反射影响是非常大的^［34］。表2的计算结果再次证明了这一点，从表中可以明显看到，各工况下鱼礁型潜堤群的反射系数均大于0.75。随着鱼礁型潜堤间距比的增加，当d/l_r=0.5~2.0，反射系数K_r随着d/l_r的增大而减小，当d/l_r=2.5时，反射系数骤增，达到峰值0.974，而当间距比继续增大时，反射系数又开始下降。在有限水深下，等间距连续布置的系列鱼礁型潜堤可以视为海底周期性起伏的地形，当表面重力波通过连续潜堤时，如果系列潜堤的间距在半倍波长的整数倍附近，波浪与地形将发生共振，波浪反射达到最强，此为流体力学中的布拉格共振现象^［35］。d/l_r在本研究中仅针对鱼礁型潜堤布设间距的测试工况，不反映潜堤与波长的关系。为进一步分析鱼礁型潜堤间距比与波浪反射程度之间的定量关系，定义相邻潜堤中心线之间的距离为S（如图1所示），即S=d+l_r。基于布拉格共振的基本思想，计算各工况下2S/L_i（L_i为波长）的数值，所得结果亦列于表2中。从表中可以看出，在间距比d/l_r=2.5工况下，2S/L_i≈1，满足布拉格反射共振的条件，说明此时波浪与鱼礁型潜堤群同频共振发生了布拉格反射，因此反射系数最大，此时透射系数几乎最小。随着鱼礁型潜堤群间距比的增大，透射系数总体上随之逐渐减小，说明透水鱼礁型潜堤的消浪效果随着布设间距的增大而有所增强。耗散系数的变化趋势与反射系数相反，但变化范围较大，当d/l_r=2.5时，K_l最小值为0.043，其余工况下耗散系数为0.330~0.416，数值上相差一个量级。本研究中给定的入射波浪属于小周期波浪，水质点运动轨迹较短，鱼礁型潜堤对波浪表层水质点的扰动也较小，在波浪与潜堤的交界处，通过潜堤透水圆孔结构进出的紊流对波浪作用明显，耗散了部分波浪能量。

表2  鱼礁型潜堤反射、透射、衰减系数

Tab.2  Reflection, transmission and dissipation coefficients of the reef-type breakwater

工况	间距比 d/lr	2S/Li	反射系数Kr	透射系数Kt	耗散系数Kl
1	0.5	0.45	0.777	0.130	0.381
2	1.0	0.60	0.776	0.105	0.387
3	1.5	0.75	0.765	0.112	0.402
4	2.0	0.90	0.757	0.104	0.416
5	2.5	1.04	0.974	0.094	0.043
6	3.0	1.19	0.813	0.090	0.330

为研究鱼礁型潜堤上波高的变化特性，基于SPH数值模型计算结果，图5给出了不同工况下波浪经鱼礁型潜堤传播过程中平均波高的沿程变化曲线。由图可知，在波浪行进到堤顶上W4（x=32.8m）测点之前，W1~W3各测点的波高由于入、反射波相位不同叠加后波动显著。当波浪到达堤顶以后，波浪发生浅化变形，波能耗散较大，波高在W4测点之后急剧减小，虽然堤顶上方自由水体的水深只有1.7cm，但由于受到透水鱼礁型潜堤内部尤其是顶板圆孔内水体以及堤顶上涌流的影响，波浪在堤顶传播的过程中始终未发生破碎，而是削弱波峰继续向前传播；波浪越过鱼礁型潜堤后，不同工况下测点W5~W7波高均趋于稳定， 沿程变化很小，各点波高基本小于0.6H₀。值得注意的是，在d/l_r=2.5工况下，堤顶波高骤升，远高于其他工况下该测点的波高值，这由两方面原因引起，主要为波浪在该工况下发生了布拉格共振反射，反射波与入射波叠加在潜堤前形成了“水墙”的结果；此外，由于堤顶水深仍处于有限水深范畴内，根据浅水变形系数的规律，在该水深范围内的波高值是增大的，但增幅很小，仅为波高的十分之一。

图5  不同工况下波高的沿程变化

Fig.5  The variation of wave heights along the reef-type breakwater under the various scenarios

3.3　鱼礁型潜堤近区波生流场变化特性

波浪与透水鱼礁型潜堤相互作用过程中引起的波生湍流对波能的耗散起着重要作用，常体现为波浪越过结构体时会出现流体分离的现象，加剧了波压和速度脉动，引起局部冲刷。对于淹没水深较大的情形（h_r>10cm），由波生流引起的波能损耗更大^［36］。由3.1节可知，当t=0.25T，即在1/4波相处，波浪在鱼礁型潜堤上发生浅水变形程度最大，堤顶上方形成的上涌流流速最大并通过堤顶透水圆孔进入结构内部，此时潜堤近区形成的波生湍流十分剧烈。为进一步了解波浪与鱼礁型潜堤群的相互作用机制，图6给出了波浪传播数个周期达到稳定波况后，当t=0.25T时，不同工况下利用SPH数值模型预测得到的鱼礁型潜堤的波生湍流场分布云图。从图中可以看出，各工况下，迎浪堤顶上均形成了高速上涌流区，当间距比d/l_r=0.5~2.0时，上涌流范围和速度随着d/l_r的增大而有所减小；当d/l_r=2.5时，上涌流流速骤升，随后又随着d/l_r的增大而减小。这与波浪反射系数的变化规律一致，说明在以间距比2.5布设鱼礁型潜堤群时，由于波浪发生了布拉格反射，入射波和反射波的叠加程度达到最大，消浪效果最强，迎浪堤上形成的波生上涌流范围和流速亦达到最大，此时对于海域内鱼礁型潜堤近区的底表层水体和营养物质的交换更具实效性。随着间距比的增大，相邻鱼礁型潜堤之间的低流速区域也逐渐扩大，流速沿程逐渐减小，波浪越过鱼礁型潜堤群到达尾端潜堤背浪面后，流速被大幅度减小，仅为堤前流速的10%左右。

图6  不同工况下鱼礁型潜堤近区的流速分布

Fig.6  Velocity contour near reef-type breakwaters under various scenarios

4 结论

本文针对不同间距布设的系列鱼礁型潜堤，分别采用水槽试验和数值模拟获得了短周期规则波在系列以周期性分布的单排鱼礁型潜堤群上传播时的波面变形、波高衰减和波生流的沿程变化。

（1）建立的SPH数值模型可以较准确地模拟鱼礁型潜堤上的波浪变形以及波高的空间分布特征。

（2）在有限水深下，等间距连续布置的系列鱼礁型潜堤可以视为海底周期性起伏的地形，结构体之间的间距对波浪反射的影响最明显，波浪透射随着间距比的增大而增强。在本文研究的鱼礁型潜堤形式下，当布设间距比为2.5时，波浪与潜堤作用产生了布拉格共振现象，反射波与入射波叠加程度最大，并在迎浪堤上形成范围和速度最大的上涌流；应用于实际工程时，鱼礁型潜堤群能够极大地减少传入近岸的波能，此时消减波浪和掩护海滩的效果达到最佳。

（3）随着鱼礁型潜堤间距比的增大，相邻鱼礁型潜堤之间的低流速区域也逐渐扩大，潜堤内的流速约为堤前流速的15%，而波浪越过潜堤后背浪面的流速仅为堤前流速的10%左右；相邻潜堤间的间距使波浪在传播过程中波峰水质点速度呈现增大—减小—小幅增加—再减小的变化趋势，连续变换的水流速度在一定程度也促进了工程海域的水体交换，这对工程海域水环境的改善具有重要意义，可为后续开展水环境评价工作提供理论参考依据。

上述研究结果为进一步研究大规模鱼礁型潜堤上的波生湍流、结构稳定性、透水孔隙率等问题奠定了基础，也可为后续鱼礁型潜堤的开发和结构设计的改进提供科学依据。此外，本文仅对波浪在同一高程布设下的单列鱼礁型潜堤的传播特性展开了讨论，实际工程中鱼礁型潜堤群往往由多列鱼礁型潜堤组成，且由于向海侧底高程实际中一般低于向陆侧，鱼礁型潜堤并非都置于同一高程上，因此波浪作用更为复杂。为提高工程的实用性和适用性，今后的研究中将对其进一步深化以探讨和完善波浪对大规模鱼礁型潜堤群响应机制的科学理论体系。

作者贡献声明

郑宇华：模型计算，数据分析，初稿撰写与修改；

匡翠萍：学术指导，论文审阅；

韩雪健：实验操作；

马悦：实验操作。