基于SPH的新型人工块体TB-CUBE设计优化和水动力特性数值模拟

彭 程，王 昊，张华庆，陈汉宝，赵 旭

(1.长沙理工大学 水利与环境工程学院 水沙科学与水灾害防治湖南省重点实验室, 湖南 长沙 410114；2.交通运输部天津水运工程科学研究院 港口水工建筑技术国家工程研究中心, 天津 300456)

人工护面块体[1]是海岸工程防浪建筑物消浪的典型措施，探究波浪作用下防浪建筑物及其护面块体的水动力特性如越浪、爬高等，具有重要的科学意义和实际工程价值。采用数值模拟方法探究这一问题具有不受时间和空间尺度限制的特点，在工程实践中日益普遍。随之产生大量波浪结构物相互作用的数值模型及算法。其中，基于光滑粒子（SРH）[2]的纯粒子数值方法，因其不受网格约束可适用于求解大变形的自由表面流动问题，如波浪破碎[3]、液体飞溅[4]、波浪与多孔介质的相互作用[5]等，而传统网格法解决上述问题均有一定困难。SРH在诞生时模拟真实海况受到计算条件的制约，但GРU并行算法的引入令计算效率大幅度提高[6]，使其开始广泛运用于海岸工程实践[7]，也使得模拟真实海况下的约1 000个不规则波的越浪场景成为可能[8]。Domínguez等[9]利用DualSРHysics与MoorDyn[10]进行耦合，模拟了规则波作用下锚定浮式结构的动力响应，其纵荡、横摇、首摇自由度与试验值吻合良好，有望适用于真实海况下的浮式海上风机、海洋石油平台、波浪能转换器等海工结构物的安全性评估。Crespo等[11]首次利用SРH探究了波浪与浮体式震荡水柱式波能转换器的相互作用情况，但因缺乏有效造波和吸波系统，在实际运用中受到限制[12]。Meringolo等[13]基于含有人工黏滞项的2维SРH模型，模拟了波浪与透空沉箱式防波堤的相互作用，模拟结果表明波压力在海墙孔处的降低，同时揭示经验公式[14]因未考虑开孔而导致波压力估算失准的问题。Monaghan[15]、Dalrymple[16]等研究了孤立波在海滩上的传播及波浪在不同糙率防波堤上的爬坡与越浪，但未考虑护面块石的作用。Altomare等[17]首次利用基于弱压缩型SРH法的DualSРHysics源代码[6,18]探究波浪与3维防波堤护面块石相互作用，模拟所得波浪爬高与试验及经验公式吻合较好。

护面块体在实际工程中种类繁多，包含传统的扭“王”字块、扭“工”字块、四脚空心方块及改进型四脚空心方块、韩国Chi Block[19]、本文TB-CUBE[20]等。各种块体在不同防波堤结构型式下的爬高和越浪差异较大，规范[21-22]所给出的经验公式均存在一定的适用范围，例如：波周期的限制，以及斜坡堤本身结构的复杂性等，工程运用有一定的局限性，需通过相应的物理模型试验和数值模拟进行确定。如：规范[22]缺少扭“王”字块体斜坡堤越浪量的多影响因素计算关系式，相关报道有限[23-24]。另一方面，为设计一种新型块体，普遍通过多组次物理模型试验确定块体最优尺寸型式及消浪特性，但物理模型试验也同样存在成本较高、效率较差、可变性低等不足；而基于广泛运用的开源代码模型DualSРHysics[25-26]搭建波浪-斜坡堤-人工块体作用的数学模型，并利用此模型进行优化设计的方法可以应对上述物理模型的不足。随着海洋结构物逐步往外海深水区发展，工程所处自然条件将更为复杂、恶劣，长周期涌浪作用于海工结构物护面块体上时可能发生冲蚀破坏[27]。而规范[21-22]中仅以设计波高作为块体重量选用的衡量标准，外海波浪多为长周期涌浪，不同于一般风浪及短周期涌浪，其波长较长，且波速较大，在水平方向上流动性很大，具有较大的能量，对人工块体稳定性的影响不容忽视。工程中往往采用依靠水体进入块体内部形成紊流效能的高透空性人工块体应对波浪的长周期问题。

为全面系统地探究新型适用于长周期涌浪条件下的高透空性人工块体TB-CUBE[20]的最优尺寸型式及其消浪特性，基于DualSРHysics模型建立波浪数值水槽，模拟规则波在铺设新型人工块体TB-CUBE斜坡式防波堤上的爬高、越浪和反射，并结合物理模型试验结果进行验证。随后，进行大量系列组次计算，分析不同细部设计尺寸下人工块体消浪的效果，在获得最优消浪效果下对应的块体设计参数的基础上，进一步探究该新型人工块体在不同波浪条件下的水动力特性，并拟合得出该新型人工块体的爬高经验公式，为该块体在实际工程中的应用提供科学依据和支撑。

1 波浪-斜坡堤-人工块体数值模型建立

1.1 TB-CUBE 新型人工块体简介

本文研究对象新型护面块体TB-CUBE由交通运输部天津水运工程科学研究院自主研发[20]。该块体以正方体为基块，包括4个1/4球形腔体、8个1/4圆柱腔体、4个1/2圆柱腔体，块体形状及各处尺寸定义示意如图1所示。图1中：R1为块体球相对半径，R2为圆柱相对半径，h为块体边长，a为块体球边长，b为块体基块边长，c为块体球体厚度，d为块体基块厚度。其消波原理主要是波浪在爬坡时从孔洞内进入块体与块体之间形成的消能室，在消能室形成漩涡、冲击和摩擦，耗散波浪能量，从而达到消波的目的。TBCUBE块体具有孔隙率高、安放适应性强、稳定性高、结构强度高等优点。

图1 新型护面块体TB-CUBE示意图Fig.1 Sketch of the new artificial block TB-CUBE

1.2 物理模型试验描述

斜坡堤新型人工块体TB-CUBE水动力特性断面物理模型试验在交通运输部天津水运工程科学研究院风浪流试验水槽中进行。模型几何比尺为1∶40，试验断面为斜坡式结构，斜坡护面为TB-CUBE块体，块体重量为156 g。下设厚度为25 mm、重量为护面块体重量1/10～1/20的块石垫层。试验模型布置如图2所示。试验中对2种水位（0.3、0.4 m）、4种波高（0.06、0.08、0.10、0.12 m）、5种周期（1.2、1.5、1.8、2.2、2.5 s）及3种斜坡坡度（1.5、2.0、2.5）的不同组合工况下TB-CUBE的爬高、越浪和反射进行了试验研究。该试验组次亦为后文数值模拟的计算波浪提供输入条件。

图2 物理模型布置图Fig.2 Sketch of physical model layout

1.3 SPH数学模型控制方程

波浪运动的质量守恒方程和动量守恒方程的光滑粒子流体动力学（SРH）离散形式如式（1）和（2）所示：

式（1）～（2）中：ρ为粒子密度；t为时间；i、j为粒子序号；M为粒子质量；V为粒子速度；P为压强；g为重力加速度；E为人工黏滞项；W为核函数，如式（3）所示：

式中，q=r/hSL，r为任意给定的两个粒子i和j之间的距离，αD= (hSL)3/π，hSL为定义光滑核函数影响区域范围的光滑长度。DualSРHysics基于弱可压缩流体假设，使用式（4）来确定流体压力：

1.4 TB-CUBE块体结构与斜坡堤数值模型建立

数学模型采用与物理模型试验相同的几何比尺，利用可视化3维AutoCAD建立TB-CUBE块体结构，以STL 文件形式导入DualSРHysics生成离散粒子。为更精细地模拟块体细部形状，经调试，粒子间距取0.001 m。斜坡堤数值模型建立方法与块体的建模方法类同，TB-CUBE 护面块体结构与斜坡堤数值模型建模实景如图3所示。

图3 护面块体与斜坡堤数值模型建模实景Fig.3 Conceptual diagram of breakwater with artificial block in the numerical model

DualSРHysics模型通过主动吸收式推波板造波程序提供动力边界条件和周期性边界条件。在3维波浪水槽模型中，由于水槽边壁与波浪传播方向一致，边壁粒子会对邻近流体粒子产生特别大的排斥力，产生不实的数值耗散，可利用周期性边界条件避免此问题，而其余边壁则采用默认的动力边界条件。3维数值波浪水槽全长8.0 m，宽0.1 m，高1.0 m；造波板宽0.1 m，高0.6 m，位于水槽左端0.5 m处。

利用构建的数值波浪水槽进行TB-CUBE块体爬高及越浪验证，以确保数值波浪水槽的可靠性。斜坡堤位于距离造波板6 m处，并在经简化后的不透水斜坡上铺设两列TB-CUBE块体，块体模型边长0.05 m，堤前布设波高仪G1和G2用于采集堤前波高；堤后布置长0.3 m、宽0.1 m、高0.4 m的集水箱采集越浪。为兼顾精确模拟块体形状与计算精度的要求，粒子间距取0.001 m，粒子数为7 678 469。模型采用规则波入射，波要素条件与物理模型试验相同，模拟时间为30 s，模型布置如图4所示。部分代表性的爬高和越浪验证结果见表1和2。

表1 数学模型和物理模型爬高试验结果对比Tab.1 Comparison of wave run-up between numerical model and physical model

图4 爬高、越浪数学模型布置图Fig.4 Layout of the numerical model for wave run-up and overtopping

由表1及表2可知，数值模拟的爬高和越浪计算值与物理模型试验值误差分别不超过6%和9%，整体较为吻合，且反射系数变化基本趋势一致。由此可见，所建立的数值模型可较为准确地模拟波浪沿TB-CUBE护面斜坡堤的爬高和越浪。

表2 数学模型和物理模型越浪试验结果对比Tab.2 Comparison of wave overtopping between numerical model and physical model

2 TB-CUBE块体水动力特性数值模拟结果与分析

2.1 块体结构优化比选分析

TB-CUBE块体是由正方体基块挖去一定尺寸的球体和圆柱体得到具有相应孔隙率的异形块体，因此，可以通过改变球体和圆柱体的尺寸，进而改变其孔隙率，利用数值波浪水槽计算其爬高和越浪，探究其水动力特性并最终找寻最优结构尺寸。为确定最优消浪效果的块体形状，在保证块体结构整体型式不变的基础上，选取具有一定代表性的6种不同细部尺寸的块体，具体参数见表3。

表3 TB-CUBE块体尺寸参数Tab.3 Size parameter of the TB-CUBE blocks

选取水深为0.3 m，波高为0.12 m，周期为1.8 s进行上述各块体的数值计算 ，将不同块体爬高、反射系数和越浪量结果分别绘于图5。由图5可知：

图5 不同块体圆柱相对半径与爬高、反射系数和越浪量关系Fig.5 Relative radius of the artificial block versus relative wave run-up, reflection coefficient and overtopping

当块体球相对半径R1相同时，爬高、反射系数随圆柱相对半径R2的增大而减小，越浪量随R2的增大而增大。以1#和2#块体为例：由于2#块体的R2大于1#块体，当块体的圆形孔径更大时，更多的水体在爬升过程中落入孔洞，波能得以快速耗散，最终表现为2#块体的爬高小于1#块体，其值分别为0.138和0.178 m；而较小的圆形孔径使得流体粒子跃出块体孔隙的“门槛”变高，更易滞留在块体空腔而不越过堤顶，最终导致1#块体的越浪量小于2#块体，其值分别为0.659和0.758 L/(m·s)。

当块体圆柱相对半径R2相同时，爬高、越浪量随球相对半径尺寸R1的增大而减小。以1#、3#和5#块体为例：随着R1的增大，块体孔隙率增大；随着孔隙率的增大，波浪在爬升过程中被消耗的能量越大，最终表现为5#块体的爬高和越浪均为三者最小，1#块体次之，3#块体最大，爬高值分别为0.168、0.178和0.188 m，越浪量分别为0.585、0.659 和0.707 L/(m·s)。

综合考虑以上块体爬高、反射系数、越浪、块体孔隙率及材料用量等因素，6#块体，即相对尺寸R1为0.475h和R2为0.315h时综合性能最优。

2.2 块体水动力特性模拟与分析

2.2.1 规则波作用下的波浪爬高

利用本文构建的数值波浪水槽，针对第2.1节块体优化比选得到的结构最优的6#块体进行TB-CUBE水动力特性模拟研究。波浪在光滑斜坡和TB-CUBE块体斜坡在t为5.25、5.45、5.75、6.15和6.45 s时刻的爬升和回落过程如图6～10所示。

图6 波浪沿堤爬升过程（t=5.25 s）Fig.6 Wave climbing along the breakwater at t=5.25 s

图6～10分别描绘出规则波作用下波浪沿斜坡堤爬高中波浪传播至堤前、流体爬升、爬至最高点波能耗尽和流体回落的完整过程。图6显示波浪第1次传至堤前，波浪呈规则形态，粒子速度分布均匀，波峰速度大于波谷速度，满足波浪理论，水底流速较小，波前面流速较小。图7显示波浪撞击坡面并开始沿堤爬升过程，此时坡面水质点流速较大。图8显示波浪沿堤逐渐爬升，流速逐渐减小，直至为0，达到最大爬高；波浪作用于TB-CUBE护面块体斜坡时，部分水体会进入护面层孔洞，波能减小，使得波浪沿堤面爬高明显小于光滑斜坡。图9显示流体在重力作用下开始沿堤回落，流速逐渐增大，此时，TB-CUBE护面斜坡对回流产生一定阻碍作用，使得其回流速度小于光滑斜坡。图10显示回落水体与下一个波浪正面相遇，在堤前破碎，破碎区流速大且紊乱，然后开始下一周期的爬高过程。

图7 波浪沿堤爬升过程（t=5.45 s）Fig.7 Wave climbing along the breakwater at t=5.45 s

图8 波浪沿堤爬升过程（t=5.75 s）Fig.8 Wave climbing along the breakwater at t=5.75 s

图9 波浪沿堤爬升过程（t=6.15 s）Fig.9 Wave climbing along the breakwater at t=6.15 s

图10 波浪沿堤爬升过程（t=6.45 s）Fig.10 Wave climbing along the breakwater at t=6.45 s

为进一步研究TB-CUBE块体的消浪特性，通过对比以上波浪在光滑斜坡和TB-CUBE护面斜坡上的爬高测定TB-CUBE糙渗系数K△。K△与斜坡护面型式有关，由块体护面斜坡和光滑面斜坡上的爬高比值得到，计算结果见表4。

表4 TB-CUBE块体糙渗系数测定结果Tab.4 Results of the TB-CUBE roughness-permeation parameters

表4中，TB-CUBE块体糙渗系数平均值为0.475，综合比较规范[21]给出的各种护面型式糙渗系数，TBCUBE糙渗系数小于0.5，优于常用护面块体型式（如块石、四脚空心方块等），表明其消浪效果较好。

2.2.2 规则波作用下的斜坡堤越浪

规则波作用下斜坡堤越浪过程大致分为波浪传至堤前、波浪破碎、波浪沿堤爬升、波浪越过堤顶、流体回落。图11为不同时刻下斜坡堤越浪模拟过程流场实景图，由蓝到红颜色深浅代表流速大小。图11（a）显示波浪在堤前破碎。图11（b）显示波浪沿斜坡面上爬，流速较大，流体爬得最高。图11（c）显示波浪爬至堤顶，流速较大，流体跃上堤面；下层流速较小，水粒子流入堤顶块体孔隙。图11（d）显示波浪越过堤顶形成越浪，但值得注意的是部分水体由于TB-CUBE结构特性产生的“门槛效应”而被滞留在堤顶的块体空腔内。图11（e）和（f）显示水体沿坡面回落，在堤前形成漩涡，流速紊乱，然后开始下一周期的越浪过程。

图11 不同时刻下斜坡堤越浪模拟过程流场实景图Fig.11 Wave overtopping and velocity distribution on breakwater at different times

3 讨 论

3.1 斜坡坡度的影响

首先分析入射波高与斜坡堤坡度对波浪爬高的影响。图12和13为相对爬高R/H（R为波浪爬高值）与坡度及入射波高的关系。由图12、13可知，R/H与坡度m呈单峰曲线趋势，即爬高随坡度变缓先增大后减小。当斜坡较陡时，大部分波能被反射消耗，此时爬高较小；当坡度m由1.5变为2.0时，斜坡变缓，反射减小，部分波能转化为势能，爬高大幅增大，以波高为4.8 m、周期为6 s为例，波浪的相对爬高由1.23增至1.67；当坡度由2.0变为2.5时，斜坡进一步变缓，波浪在堤前破碎，耗散大量波能，此外，波浪在缓坡上爬升距离加大，沿程摩阻和护面层中水流紊动的波能耗散作用增强，因此爬高变小，同样以波高为4.8 m、周期为6 s为例，波浪的相对爬高由1.67降至0.99。

图12 不同坡度和波高、周期条件下的相对爬高（d=16 m）Fig.12 Slope versus relative wave run-up under various wave heights, periods at d=16 m

图13 不同坡度和波高、周期条件下的相对爬高（d=12 m）Fig.13 Slope versus relative wave run-up under various wave heights, periods at d=12 m

3.2 波坦的影响

图14为不同坡度下反射系数、相对爬高与波坦的关系。由图14（a）可见，反射系数随斜坡堤坡度的变缓而减小。当坡度一定，波坦（波长/波高）较小时，破碎波浪占主导；随着波坦增大，波浪破碎减小，波浪反射逐渐增大，当波坦达到一定值时，反射系数趋于稳定。图14（b）和（c）为不同周期时相对爬高与波坦的关系图，由图14（b）和（c）可见，相对爬高（R/H）与波坦（L/H）总体呈正相关，即相对爬高随着波坦的增大而增大。波坦较小时，大量波浪在堤前破碎，波能消耗；随着波坦逐渐增大，波浪破碎程度随之减小，耗能减小，因此前期爬高值急速增长；但随着波坦进一步增大，破碎波进一步减少，此时波浪反射加强，反射波能增加，爬高值呈缓慢增长。

图14 波坦与反射系数、相对爬高关系Fig.14 L/H versus reflection coefficient and wave run-up

3.3 波浪破碎的影响

采用广泛运用的破波参数Iribarren 数ζ0[28]来描述波浪与斜坡堤相互作用时的波浪破碎形态，以判定破碎波形式属于崩破波（ζ0≤0.5）、卷破波（0.5＜ζ0＜3.3）或激破波（ζ0≥3.3）。图15为相对爬高与破波参数之间的关系。由图15可知：当m=1.5时，最小破波参数2.33，波浪多为激破波；当m=2.0，且当波高大于3.2 m时，波浪表现为卷破波，波高小于3.2 m时，主要为激破波；当m=2.5时，破波参数在1.40～3.31范围内，破碎波浪均属于卷破波。总体而言：相对爬高随破波参数的增大而增大，当破波参数小于2.5时，此时波浪主要表现为卷破波，即波浪完全破碎，反射波较少，爬高与破波参数表现为良好的线性关系；当破波参数大于2.5时，相对爬高随破波参数的增长变缓，这是因为，随着破波参数变大，波浪在堤前以激破波的形式破碎，再加之与反射波叠加效应，导致波浪爬高与破波参数呈非线性增长。

3.4 爬高公式拟合

以无量纲参数R/H和相对水深L/d为参数，对不同坡度下的R/H与L/H之间的关系进行非线性拟合，以此来确定斜坡堤不同坡度下的对应的参数，拟合形式见式（5），其中，波浪爬高拟合公式系数A(m)、B(m)、C(m)的值见表5。

表5 波浪爬高拟合公式系数Tab.5 Fitting coefficients for the wave run-up formula

进一步把表5中不同坡度下的系数A(m)、B(m)、C(m)与坡度m值进行拟合，分析得到坡度函数与坡度之间的关系式。通过将各个公式进行拟合后，最终得到波浪在新型人工块体TB-CUBE上的爬高经验公式，如下所示：

为了分析式（5）的精确度，将数学模型试验值与公式计算值进行对比，如图16所示，采用最小二乘法进行回归分析，得到相关系数为0.981。由图16可知拟合公式计算值与试验值线性相关性很强，说明本文拟合爬高公式具有较高精确度，能够真实反映爬高情况。

图16 数学模型试验值与拟合公式计算值对比Fig.16 Wave run-up from model results versus fitting formular

4 结 论

1）基于DualSРHysics代码可以模拟波浪在铺设新型高透空性护面块体TB-CUBE防波堤上的水动力特性。优化比选结果表明：当块体球体半径R1一定时，爬高、反射系数随圆柱尺寸R2增大而减小，越浪量随圆柱尺寸R2的增大而增大；块体圆柱半径R2一定时，爬高、越浪量随球体半径R1的增大而减小。

2）对优化后的TB-CUBE块体展开水动力特性数值模拟研究，所得糙渗系数平均值为0.475，表明该块体消浪效果良好。TB-CUBE块体本身的结构特性对波浪与斜坡堤的相互作用过程影响很大，护面层圆孔对波浪的爬高和回落均有阻滞作用，块体内部孔隙可以容纳水体，形成紊流效能，其“门槛效应”也可以相应地减少越浪。

3）波浪沿TB-CUBE斜坡堤爬高数值模拟结果表明，波浪爬高随波高、波坦增大而增大，随坡度的变化呈现单峰曲线趋势。并进一步拟合了爬高公式，可供相应工程参考。