基于SPH 方法的沉箱开孔位置对消浪效果的影响

刘枫，甘忠勇，王坤，韩丽影，谷雪清，金瑞佳

（1.长江南京航道工程局，江苏 南京 210000；2.河海大学港口海岸与近海工程学院，江苏 南京 210098；3.交通运输部天津水运工程科学研究院，天津 300456）

0 引言

我国是海洋大国，随着海洋工程领域建设的不断发展，传统实体直墙式的防波堤结构已经逐渐不能满足实际工程的各种需要。近年来具有低反射、低越浪、受力小和工程造价低等优点的开孔沉箱结构，在海洋工程的实际建设中得到广泛应用。

目前对于开孔沉箱与波浪作用的研究主要是通过物理模型试验和数值模拟。其中物理模型试验因其得到结果比较可靠、直观，成为开展开孔沉箱研究的重要手段之一。

传统的开孔沉箱物理模型试验大多基于弗劳德相似准则进行[1-5]。但由于前壁孔洞的存在，开孔沉箱受流体黏性的影响显著，有学者研究发现按照弗劳德相似准则设计的开孔沉箱的试验研究存在比尺效应问题[6]。而后有学者采用物理模型对开孔沉箱模型的比尺效应进行了研究[7]，发现了开孔沉箱的内部受力与入射波要素的关系与模型比尺相关。然而物模试验的研究依赖于不同比尺的试验设施，成本高昂，因此很多学者开始采用数值模拟的方法研究开孔沉箱结构的相关问题。如Liu 等通过建立数值水槽，分析了波浪与开孔沉箱的相互作用[8]、波浪的反射和透射[9]以及斜向的波浪影响[10]。

随着SPH 计算方法的发展，一些学者开始验证并使用SPH 方法研究开孔沉箱相关问题[11-12]。如唐晓成[13]、任喜峰[14]和龚凯[15]等，其数值结果与物理模型试验结果吻合良好，证明SPH 方法对开孔沉箱结构研究的可行性。SPH 方法是一种拉格朗日形式的无网格数值计算方法，适用于模拟黏性流体强非线性运动。在波浪作用下，开孔沉箱孔洞附近自由液面复杂、湍流运动剧烈，SPH 方法的优势得以很好体现。因此本文采用SPH 方法对开孔沉箱相关问题进行研究。

而关于开孔沉箱的研究成果主要集中在其消浪效果上。Tanimoto 等[16]采用物理模型对开孔沉箱进行消浪效果的分析和研究，指出了当沉箱消浪室宽度为波长的0.18 倍时，其消浪效果最为明显。Suh 等[17]进行了开孔沉箱的物理模型试验，提出了当消浪室宽度为入射波长的0.2 倍时，反射系数Kr最小。陈雪峰等[18]通过研究总结了反射系数Kr与消浪室宽度和开孔率的关系。万庆宇[19]、路伟[20]分别在规则波和不规则波的作用下，通过改变沉箱基床高度，拟合了反射系数Kr与开孔率、消浪室相对宽度、相对水深等影响因素之间的关系式。

通过上述研究成果可以发现，消浪效果主要与开孔率、消浪室相对宽度、相对水深d/L 等因素相关，然而关于不同开孔位置对开孔沉箱消浪效果的影响研究成果还较少，而沉箱开孔位置也是影响其消浪效果的重要因素之一，所以本文将采用SPH 方法，对不同开孔位置的开孔沉箱的消浪效果进行研究，分析不同开孔位置对开孔沉箱消浪效果的影响，并提出使开孔沉箱的消浪效果达到最佳的开孔位置设计方式，对开孔沉箱实际工程中消浪研究具有一定参考意义，所得成果可为相关未来工程建设提供借鉴。

1 基本理论

1.1 SPH 形式的控制方程

流体运动方程由连续性方程和Navier-Stokes方程控制，在弱可压缩SPH 模式[21]中的连续性方程和Navier-Stokes 方程为：

式中：下标“i”和“j”分别表示目标和相邻的粒子；ρ 为密度；V 为粒子体积；u 和r 分别为速度和位置向量；表示i 方向的空间向量；ɡ 为重力加速度；Wij为W（ri-rj，h）核函数的简写形式[22]；h为光滑长度，h=1.5dp；dp为粒子间距；δ 为取为0.1 的系数；c0为数值声速[23]；ψij为用来稳定数值压力场的耗散项[21，24]；α=8v/（hc0）为人工黏性系数[25]；v 为运动黏度；ρ0=1 000 kg/m3为参考密度；πij表示速度向量和位置向量的相对值；ψij和πij表达式分别如下：

式中：⊗表示张量积，若u 的导数从式（2）中获得，则流体粒子的运动位置可以通过式（6）得到：

1.2 模型参数设置

1）造波方法

本次数值波浪水槽的左侧造波区采用主动吸收式造波方法，可以降低二次反射波浪对入射波的影响。采用的主动吸收式造波技术为Hirakuchi等[27]的时域控制及方法。

2）边界条件

本模型的固体壁面边界由改进的动力边界粒子法[28]施加。在固体结构的外轮廓上布置两排动力边界粒子（简称DBPs），DBPs 密度用连续性方程（1）计算，同时计算DBPs 支持域内流体粒子的平均密度，取两部分计算结果的加权平均值作为DBPs 的修正密度，修正公式如下：

3）时间步长

采用可变时间步长，其计算公式[29]如下：

式中：等号右侧第一项是对粒子运动加速度的限制，系数ξ1=0.25；等号右侧第二项是对流体黏性扩散的约束，系数ξ2=0.125；等号右侧第三项是融合流体黏性的CFL 条件，系数ξ3=0.5。

对于反射系数的计算，本文采用两点法对入射波、反射波进行分离。利用两点浪高仪处的波面高度变化信息，通过希尔伯特变换得到波浪信号在复域内的解析表达式，从而达到波浪信号的分离。根据波浪信号分离的结果，得到入射、反射波的幅值以及反射系数等参数，进而得到它们之间的相位信息。

2 数值水槽模型建立

基于SPH 方法建立数值波浪水槽，其计算域如图1 所示。水槽上游布置主动吸收式推板造波机[30]，可以避免沉箱的反射从而生成稳定的波浪，水槽下游布置开孔沉箱。造波板距离原点0.2L，水槽长度为2L，L 为试验波长，开孔沉箱的宽度为B，水深为d，波浪水槽由左侧造波区、底板和右侧开孔沉箱3 个固体壁面边界组成。采用开孔沉箱纵向断面代表开口沉箱模型，孔形为横条孔，沉箱开孔高度分别为a1和a2，模型开孔率μ=，即为开孔总高度除以沉箱高度，沉箱外取3 个浪高仪WP1—WP3，间隔为0.1L。

图1 波浪作用于开孔沉箱问题计算域示意图Fig.1 Schematic diagram of calculation domain for problem of wave action on perforated caissons

2.1 模型验证

为了验证模型的准确，将应用本模型计算得到的波浪反射系数与Liu 等[6]的解析结果进行比较，以验证此模型的可靠性。取计算域水深d=0.22 m，入射波高H=0.03 m，波周期T=1.0 s，开孔率μ 取20%和30%两种。由图2 可以看到，模型计算与学者解析得到的波浪反射系数吻合较好，且两者随舱室相对宽度B/L 的变化趋势一致，从而证明可采用本文建立的SPH 模型进行相关分析研究。

图2 本模型的反射系数计算结果与解析解的对比Fig.2 Comparison of reflection coefficients calculation results of the modal and analytical solutions

2.2 模型计算设置

模型计算示意图同图1，开孔沉箱宽度B=0.13L，水深h=0.5 m，计算过程中保持L3不变，通过改变上部开孔位置即通过缩短L1同时加长L2进行相关研究，模型开孔高度a1=a2=0.08 m，得到开孔率μ=20%，试验中开孔率μ 保持不变，沉箱上部开孔的顶端在水下0.08 m 称为开孔位置1，上部开孔的顶端在水下0.04 m 称为开孔位置2，开孔的顶端与水平面齐平称为开孔位置3，如图3所示，在4 种波陡kA 下开展3 个波浪周期入射波的数值计算，具体工况参数如表1 所示，表中H 为入射波高。对比不同开孔位置的开孔沉箱的反射系数Kr、波面高度变化、涡量场和流场变化分析其消浪效果。

表1 开孔沉箱算例的计算工况参数Table 1 Calculation working condition parameters of perforated caisson examples

图3 不同算例模型设置图Fig.3 Sketch of different calculation models

3 开孔位置对消浪效果的影响

3.1 不同开孔位置的反射系数对比

反射系数Kr是反映消浪效果的一个重要参数，所以对比不同开孔位置的反射系数Kr来反映消浪效果，不同位置的反射系数Kr随kA 变化如图4 所示。可以看出，各个周期波浪入射下不同位置的反射系数规律相近，对比3 种开孔位置下的反射系数发现，第1 种开孔位置反射系数Kr最小，即消浪效果最好；第2 种、第3 种开孔位置的反射系数较接近，均大于第1 种开孔位置。

图4 反射系数随开孔位置的变化规律Fig.4 Variation law of reflection coefficient with different open hole location

3.2 不同开孔位置的波面高度变化对比

开孔沉箱与传统直立沉箱相比，波浪可以通过开孔进入舱室，不同开孔位置对波浪通过开孔进出舱室的作用方式有影响，从而影响沉箱外的波面高度变化。沉箱外3 个浪高仪与造波板的距离分别为1.1L、1.2L 和1.3L，由于本文算例较多，仅展示波浪周期T=1.2 s，波高H=0.03 m 时的波面高度η 时间历程曲线，如图5 所示。

图5 不同开孔位置的波面高度变化历时曲线Fig.5 Duration curve of wave surface elevation variation with different open hole locations

从图5 看出第1 种开孔位置波面高度变化最小，消浪效果最好，第2 种、第3 种开孔位置的波面高度变化曲线呈现无规律。第2 种开孔位置的消浪效果比第3 种开孔位置好，即第2 种开孔位置波面高度变化小于第3 种开孔位置波面高度变化，如图5（a）和图5（b）；第2 种、第3 种开孔位置消浪效果相似，即2 种开孔位置的波面高度变化相似，如图5（c）。

3.3 不同开孔位置的涡量场对比

在波浪与开孔沉箱相互作用过程中，孔洞周围的流体黏性效应显著，开孔沉箱受到流体黏性作用，更多波能以湍流形式耗散，进而影响其消浪效果。入射及反射在开孔处的边界层内黏性力消耗的能量越多，产生的涡量越大，消浪效果越好，其中涡量值由式（9）计算：

式中：m 为水质点质量，图6 为不同开孔位置一个周期的涡量场对比，观察不同开孔位置的开孔沉箱开孔处涡量消耗的大小，反映出不同开孔位置对消浪结果的影响。展示对比算例的入射波周期T=1.1 s、入射波高H=0.03 m，选取沉箱开孔附近的流场进行分析。

图6 不同开孔位置一个周期的涡量场对比图Fig.6 Comparison of vorticity fields in one wave period with different open hole locations

由图6 可看出，涡量场的差异主要集中在上部开孔位置，第1 种开孔位置进行入射（t=0.25T）和反射（t=0.75T）都消耗大量的能量，涡量消耗的总量最大。不同时刻的3 种开孔位置的涡量场最值不同，与入射和反射运动相关。

当波浪入射时，沉箱外水粒子流入沉箱消耗能量，第1 种开孔位置涡量的最值为27.88 m2·s，第2 种开孔位置涡量的最值为19.83 m2·s，第3 种开孔位置的最值为14.90 m2·s，第1 种开孔位置的涡量场能量消耗最大；当波浪反射时，沉箱内水粒子流出沉箱消耗能量，第1 种开孔位置涡量的最值为-30.85 m2·s，第2 种开孔位置涡量的最值为-23.25 m2·s，第3 种开孔位置的最值为-19.37 m2·s，同样第1 种开孔位置的涡量场能量消耗最大。因为第1 种开孔位置在水下的开孔上壁比较长，入射及反射在开孔处的边界层内黏性力消耗的能量较多，产生的涡量较大，消浪效果较好。

3.4 开孔沉箱附近的流场对比

与传统沉箱相比较，开孔沉箱降低反射的效果与波浪出入消浪室过程密切相关，流场是反映消浪效果的重要参数，消浪效果越好，沉箱开孔处的流速变化越剧烈，沉箱内的流速越平稳。下面对比不同开孔位置在入射波周期T=1.1 s、入射波高H=0.03 m 波浪条件下的速度场V，对比图如图7 所示。

图7 不同开孔位置一个周期的速度场对比图Fig.7 Comparisons of flow velocity field in one wave period with different open hole locations

通过图7 可以看出在开孔处及舱室内，第1种开孔位置的流场变化最大，在进行波浪入射（t=0.25T）和波浪反射（t=0.75T）时，3 种开孔位置流场的流速较大，第1 种开孔位置的流场变化值是其他2 种开孔位置的1～2 倍，同时沉箱内流场相对最为平缓，说明大量流体在沉箱内速度发生了相互抵消，进一步证明了第1 种开孔位置的消浪效果较好。

4 结语

本文基于SPH 方法建立了波浪与开孔沉箱相互作用的数学模型，保持开孔沉箱下部开孔的位置不变，只改变上部开孔的位置，开展了不同波浪要素下数值计算分析。通过对比不同开孔位置的反射系数Kr、波面高度变化、涡量场和流场变化，得到以下结论：

1）当保证波浪能充分进入开孔沉箱内部情况下，沉箱上部开孔的顶端距离水面最大时，反射系数最小，消浪效果最好，随着开孔位置的提高，反射系数和波面高度变化呈现无规律性且较大。

2）当沉箱上部开孔的顶端距离水面最大时，在沉箱开孔处的涡量消耗的总量最大，消浪效果最好，且在t=0.5T 时，沉箱前墙前后侧波面存在较明显的水位差，说明此种开孔位置的开孔沉箱消浪效果明显。

3）当沉箱上部开孔的顶端距离水面最大时，其开孔处流速变化值最大，是其他2 种开孔位置的1～2 倍，沉箱内流速最为平稳，进一步说明此种开孔位置的沉箱的开孔处及舱室内流场变化剧烈，消浪效果最好。