木本植被覆盖岸坡上波浪爬升过程的数值模拟研究

徐海珏，胡萍，白玉川*，杨波

(1.天津大学 水利工程仿真与安全国家重点实验室，天津 300350；2.天津大学 河流海岸工程泥沙研究所，天津 300350；3.AECOM(天津)工程顾问有限公司，天津 300072)

1 引言

近年来，海平面上升、海岸线侵蚀等问题对于沿海地区产生了不利影响[1]。伴随着沿海住宅和生态景观的建设开发，人们愈发强调海岸防护措施的功能性和生态友好性。近岸木本植物构成的生态缓冲带作为一种新型的生态防护结构，不仅能建立起近岸的绿色生态系统持续改善环境，还可以显著地消减波能，减少波浪对近岸结构物的影响，更能与传统的近海工程防护措施共同形成工程植物防护体系[2-3]。基于这些优势，利用木本植物形成的植被缓冲区的生态软防护结构正在逐渐取代传统的硬质海岸保护结构。因此，分析植被和波浪因素对消浪效果的影响，深刻理解波浪－植物之间相互作用的机理，如何准确预测波浪在植被影响下的传播过程，提升其对海岸的防护效果是当前亟待解决的问题。

近年来，针对植被消浪的研究日渐深入，且学者多采用理论推导、试验分析和现场调研3种方式。相对于理论推导，物理缩比试验和现场观测研究可给出针对具体工程的波浪变形现象的定性描述。在物理模型试验方面，白玉川等[4]选取桧柏树枝作为防浪树的试验模型，探究了防浪林种植宽度、底坡高度、底坡高程对消浪效果的影响；吉红香[5]用刚性树干和带有柔性枝叶的模型树模拟堤外滩地防浪林，分析了滩地水深、浪高、波长等因素对植物消浪的影响，以及波高、波浪爬高等的分布规律；周悦等[6]进行刚柔组合型植被消浪特性的物理试验研究，探究了防浪林的树木大小、刚柔等不同组合方式对消浪效果的影响；Yao等[7]通过PVC直管模拟植被冠层，研究推导刚性植被阻力系数的方法。在现场观测方面，Mazda等[8]和Quartel等[9]通过对波浪入射条件和传播过程中的波浪变形的分析，分别讨论了越南Tong King三角洲和红河三角洲造林对海岸保护的作用；曹大正等[10]、李志强等[11]和田野[12]基于现场的波浪和泥沙的动态观测，分别研究了互米花草、红树林、无瓣海桑等近岸植物的消浪促淤效果。尽管有很多学者做出了创新的成果，模型试验和实地调研仍有其各自的局限性：实验室模型相似性选择的确定的合理性有待进一步讨论，且受到很多实验室客观条件（动力、侧限、波浪反射）的影响；实地调研需要大量的人力物力，还会对环境产生一定程度的影响。由于波浪爬坡过程是一较为复杂的非线性过程，其演变过程受多方面因素的影响，因此有必要通过理论推导和数值模拟对其进行讨论，并给出沿程波面和流场变化的细致的描述。在理论推导方面，Bai等[13-14]、李绍武等[15]从N-S方程、大涡模拟等理论出发，建立了波浪浅水传播变形模型；Dalrymple等[16]基于线性波和刚性圆柱体空间分布的摩擦模型对局部波浪耗散进行了计算；Fernando和Losada[17]基于理论推导和经验系数相结合的方式，用刚性直杆代替植被，给出了可变深度植被场作用下波浪传播模型；杨建民[18]引入反映刚性植物树干拖曳力作用的水流阻力系数λ，得到波动势函数的一阶和二阶近似解；唐军等[19]基于考虑波浪折射、绕射等综合效应的抛物型缓坡方程，建立了近岸植被影响下波浪传播运动的数学模型；Maza等[20]基于IHFOAM的三维数值方法来研究海啸波与红树林的相互作用，并用拖曳系数的宏观方法计算波浪力；Zhu和Chen[21]讨论了由植被阻力引起的耗散和近共振三位一体的相互作用对单个谐波衰减的影响。但目前的理论描述多将植被简化为刚性直杆，缺乏考虑植被枝干的分区作用和植被密度、树干高度和树枝角度等植被特性对波浪的影响，且对在倾斜底坡这一条件下的植被的消浪效果缺乏描述和分析。因此，本文参考了以海南岛文昌海岸近岸的植被（图1）为代表的近岸木本植物形态，建立了木本植被保护下波浪沿斜坡爬升的表面波衰减理论模型。针对植被上、下层对波浪不同的作用形态，分别引入相应的阻力系数以反映植被树枝和树干对波浪的拖曳力作用，并采用有限差分法离散方程和标记网络法（Marker and Cell Method，MAC）来跟踪自由曲面上的水颗粒轨迹。以波浪沿坡度为1/30的斜坡爬升为算例，分别讨论了有无植被作用下沿程波形和流场的变化，并分析了植物枝干的直径、高度、密度、角度等植被特性和波浪因素对消浪效果的影响，得到植被消浪的基本规律。

2 理论模型

理论模型将波浪简化为垂直平面上的二维波，得到波浪的连续性方程和修正的Navier-Stokes（N-S）方程，并引入阻力分量来反映植被枝干对波浪的作用。

2.1 控制方程

参照图2所示的坐标系，水平x坐标以向右作为其正方向，垂直z坐标以向上作为其正方向，以i=tanθ为底坡斜率，θ为坡角。考虑单向波通过木本植被缓冲区传播到岸边的情况下，波浪会同时受到植被和倾斜底床的作用，木本植被海岸示意图如图3所示。

在实际的海岸通常有一个从深海到海岸的斜坡，数值模型的分析区域包含了一段平坦的河床和一段倾斜的河床。以线性减小水深为例，水深公式为

模型在水平方向上可分为3个区域讨论。区域I（-∞＜x≤x0，-d0≤z≤0）是在有植物作用之前的平坦河床部分，这个区域的入射波在平底上传播且区域中的反射波主要来自植物反射。区域II（x0＜x≤x1，-d(x)≤z≤0）则被有着不同密度和直径的树干的木本植物所覆盖。波浪从区域I传播通过植物到达区域III，而有一些反射波从区域III开始传播通过植物到达区域I，由于植物树干与波浪之间的相互作用，波能被大量耗散，即波浪衰减。区域 III（x1＜x≤xb，-d(x)≤z≤0）是在植物之后的区域，入射波的动能已经大大减弱。有时，流速降至远低于沉积物起动速度的值，以致于大多数泥沙颗粒无法沿着斜坡移动。

图1 海南岛文昌海岸近岸木本植被缓冲带Fig.1 Nearshore woody vegetation at Wenchang coast in Hainan Island

图2 无植被海岸示意图Fig.2 Sketch of beach without vegetation

图3 木本植被护岸示意图Fig.3 Sketch of beach with woody vegetation

一般来说，在没有植被的情况下，破波点存在于在区域II或区域III中的某个地方。但由于植被的存在，区域II的消波率发生了很大的变化，引起流场的变化，并最终导致了破波点位置变化。在这3个区域上，波浪可简化为二维波并提出了考虑波浪传播整体效果的连续介质等效模型，得到对应的连续性方程和修正的N-S方程：

式中，u、w分别为x、z方向的速度分量；t为时间；ρ为水的密度；p为压力；ν为运动黏滞系数；νt为动态涡流黏度；g为重力加速度；Fx和Fz分别是植物树干对波浪作用力在x方向和z方向的分力。

Massel等[22]已经证明了植物树干对流体的作用力主要为阻力，因此，Fx和Fz可以直接写成阻力分量，而Fx和Fz在不同区域的表达是不同的，具体如下：区域I：如果流体满足非黏性、不可压缩假定，则表达式为

区域II：能量耗散主要是由于流体质点与底部的相互作用，以及在没有植物的倾斜底坡上攀爬的波浪的流体质点之间的相互作用。如果在流场中增加了木本植被作用，那么枝干与流体之间的相互作用也就增加了，具体表现为枝干的阻力。由于前两种情况已经在许多研究中被频繁地讨论，本研究主要聚焦于最后一种情况，即增加的枝干的阻力。此时，仍假设流体满足非黏性、不可压缩假定，但植被的阻力不可忽略，则表达式为

首先，考虑单个小区域的木本植物作用，并将单位面积植株的数量记为N。然而，即使对同一种树，枝干的密度和直径在上层或下层是不同的。一般的木本植物上层有着浓密且直径较小的树枝，而下层为直径较大的树干。因此，我们在该模型中分区引入了不同的阻力表达式来区分树干和树枝的作用。

在植物下层 [-d(x)≤z＜-(d(x)-d1)]，单位面积上树干的数量NL明显等于植株的数量N，记d1为对应植物平均树干高度（hL）的下层水深。假定第i个植株的树干直径为DLi，那么拖曳力表示为

式中，Flx和Flz分别是下层植物树干对波浪作用力在x方向和z方向的分力；Cdmi是第i棵植株的拖曳力系数；DLi是第i棵植株的树干的的直径。

在植物上层[-(d(x)-d1)≤z≤0]，单位面积上树枝的数量NU要大于植株的数量N。此外，并不是所有的树枝都是垂直于静水面的，其阻力表达式定义为

式中，Fux和Fuz分别是上层植物树枝对波浪作用力在x方向和z方向的分力；θi是树枝与垂直方向的夹角大小，一般满足正态分布；DUi是第i棵植株的树枝的直径。

根据Massel等[22]的假设，将这些角度均用平均角度θm来表示时，理论结果与实验结果相差不大。因此，为简化模型，用平均角度θm代替各个树枝与垂直方向的夹角大小θi，得到简化后的植物上层阻力表达式定义，

如上文所述，Cdmi是一个修正系数，它不仅包括树干和树枝对流体的阻力，还包括不同圆柱形树干和树枝后面的涡旋之间的相互作用力。假设植被的生长相对较为稀疏或树干和树枝对流体的阻力起主导作用，则可以在计算中忽略涡旋之间的相互作用，那么对单个树枝而言系数Cdi与雷诺数Re相关。根据SPM[23]得到的结果，Cdi的表达式可以写成

式中，Re为雷诺数；Cdi是单棵植株的拖曳力系数。

引入以下影响因素：

式中，I(Re,N)是由于植株密度改变引起的影响因素，它的值是由另一个小范围数值模型的结果决定的。

区域III：这一区域内没有树的作用，即不存在流体和树枝之间的相互作用，且经过植被的消能作用后，在此区域的波浪形态表现为前进波和破碎波。因此，在区域III中流体不应被视为非黏性的，而应视为黏性的旋转湍流，

式中，νt是动态涡流黏度，νt=2C(x,z,t)Δ2|S|。其中，|S|为应变率的绝对值，其公式表示为|S|=[2(Sij·Sij)]1/2，Sij=[(∂ui/∂xj)+(∂uj/∂xi)]/2，Sij可表示为Sxx，Sxz，Szx，Szz，i=x,z且j=x,z；Δ为与G相关的特征过滤值；C为与湍流结构密切相关的参数；G为网格过滤值。

2.2 边界条件

自由表面边界介于剪切层和无剪切环流间，通常假设外部流完全不受湍流影响，因此在这个边界上，所有的湍流应力和通量都可以设为0。底边界按可滑移条件处理，且满足垂直于边界的速度为0和边界网格内的速度散度为0两个条件。底坡与水面相交点附近的网格既含有底边界又含有自由表面其处理稍为复杂，但因近岸岸滩和本文算例的底坡均很缓，则含有底边界和自由面的网格的压力可近似等于静水压力，以简化这些网格的处理。

3 数值求解方法

本文的数值求解方法采用有限差分法离散方程，并用MAC法来追踪自由曲面上的水颗粒轨迹。每个物理变量的网格划分和分布如图4，其中dx、dy分别为x和z方向的空间步长。

图4 MAC法单元图Fig.4 MAC method element diagram

3.1 数值方法求解流场

MAC方法采用了网格交错技术，将压力、湍流剪切应力等变量定义于网格中心，将速度 (u,w)等变量定义于网格边线中点，这样便于满足单元内速度散度为0这一条件。在方程中，对时间采用向前差分，对空间采用中心差分。由动量方程导出的差分方程以及由连续性方程导出的差分方程即为MAC方法需求解的差分方程。

求解泊松方程采用松弛迭代法，其迭代关系式为

式中，上标m和m+1分别表示第m和m+1次迭代；ω是松弛因子

其中，

式中，上标m+1/2表示由m和m+1两个迭代步导出的量。

3.2 施主格式的稳定性分析

施主格式数值耗散比中心差分格式更小，为了满足稳定性条件以及减少差分格式的数值耗散，本研究一律采用施主格式。为得到方程的收敛解，需对NS方程导出的差分方程和压力泊松方程的差分方程的数值稳定性进行分析。N-S方程为非线性方程，其稳定性分析是相当困难的。根据Von Neumann法分析，对于波浪运动，直角坐标系下二维问题的MAC方法的时间稳定性条件可简化为

泊松方程的系数矩阵是对称正定矩阵时，取松弛因子 0 ＜ω＜2，松弛迭代是收敛的，且存在着使收敛速度最快的最优 ω*值，一般来说，取 ω ＞1的超松弛迭代比取 ω ＜1的欠松弛方程迭代收敛速度要快。但在NS方程和压力泊松方程的耦合问题中，每一时间步泊松方程系数矩阵是不相同的，因此无法分析其收敛条件，也无法确定最优的 ω*值。在实际计算中欠松弛迭代法较为适用，即求解压力泊松方程时取前一时间步长的压力值作为迭代的初始值，每一步时间增量很小，压力变化也不大，因此可以保证泊松方程的迭代过程稳定地趋于收敛。

4 数值求解算例及消浪效果影响规律分析

数值求解是根据给定的边界条件，利用MAC方法对波浪的传播和爬升过程进行模拟。由于模型存在一定的局限性(恒定水深、常浪入射条件、植物均匀分布)，本文暂不进行理论模型结果与实际木本植被缓冲带观测值的比较，这二者比较的结果将随着研究的深入进行后续讨论。但为给相关的物理模型试验研究提供对照验证，也为后续物理试验的开展提供依据，本节将给出数值求解模型的算例结果，并对木本植物群的消浪效果影响规律进行分析。

本文以波浪沿1/30的斜坡爬升为算例，通过给出入射波的基本参数（入射水深DS，入射波高HS，入射波长LS，入射波周期TS），分别讨论有无植被作用下的波形变化和速度场情况。此外，在考虑木本植被作用的情况下，通过引入不同的阻力系数，讨论了植物枝干的直径、高度、密度等因素对波浪能量衰减的影响。此外，本节将得到的部分算例结果与相关文章中的试验结果规律进行对照，验证了数值模型的有效性，可供实际参考。

4.1 有无植被作用下波浪沿倾斜海滩传播的数值模拟研究

波浪沿斜堤或斜坡向上传播时，不断向上爬升必然会出现波形的不断变化，本小节通过对比分析有无植被作用情况下的波浪变形和流场情况，讨论近岸植被对波浪衰减的影响。

4.1.1 入射波浪条件

数值计算区域参照图2所示，x=0 m为波浪入射边界，给定计算水平长度xb=32 m，其中平坦河床段x0=2.0 m，倾斜河床段xb-x0=30 m，i=1/30。在初始时刻计算区域内的水体是静止的，即初始速度场为0 m/s，初始压力场按静水压力场给出。以二阶Stokes波为例，给定的各项波浪参数为：入射水深DS=1.0 m，入射波高HS=0.25 m，入射波周期TS=2.9 s，入射波长LS=8.5 m和波陡 δ =0.029。波浪通过入流边界处的水质点速度和波面位置按相对应的理论公式变化而产生。根据时间稳定性条件，网格尺度取dx=0.125 m，dy=0.025 m，时间步长dt=0.01 s，网格的划分和时间步长的选取综合考虑了差分格式的稳定性、尽可能小的数值耗散和合理的计算工作量。

4.1.2 无植被作用情况下波浪传播数值模拟及模型合理性验证

通过调整入射波波前与静水面相交点到波峰的距离，进而改变入射边界的波面位置，得到6条波浪在1/30的底坡上传播的波形曲线（图5），并以左边界位置的波面高度HL=0.03 m为例绘其传播流场（图 6）。

通过波浪传播的波形图可以发现，在没有植被的拖曳力作用时，破浪爬坡过程波高上升趋势明显，随着水深变浅，波高增大，波长变短，波前峰变陡，而波峰后逐渐出现波面变化较小的台阶式爬升、继而破碎现象，这种水平不对称的变化与流函数波理论解是一致的。对应的，以左边界位置的波面高度HL=0.03 m为例生成的流场图中，随着波前峰趋于变陡波峰处的水平速度不断增大，竖向速度大小变化加剧，并在最大水平速度达到波速时发生破碎，波浪沿斜坡传播流场（图6）不同位置处的局部放大图如图7所示。

4.1.3 有植被作用情况下波浪传播数值模拟

图5 无植被作用下波浪沿斜坡（i=1/30）传播波形Fig.5 The wave form of the wave propagation along slope (i=1/30) without vegetation

图6 无植被作用下波浪沿斜坡（i=1/30）传播流场Fig.6 The flow field of the wave propagation along slope (i=1/30) without vegetation

图7 无植被作用下流场局部放大图（流场所处的位置见图6）Fig.7 The partial enlarged drawing of the flow field without vegetation (the positon of the flow field see Fig.6)

考虑植被作用情况下，参数Cdmi设定为：Cdmi=0.7，为与无植被作用情况进行初步对比，将单位长度内的植被参数设为：树干数NL=1，树枝数NU=6，树干高度hL=0.6 m，树干直径DL=0.2 m，树枝直径DU=0.04 m，树枝倾斜度 θ=π/6，植被水平方向覆盖宽度lT=27 m，得到波浪在植被均匀覆盖的1/30的底坡上传播的波形（图10）和流场（图11）。其中，图11为布置了均匀分布的植被波浪传播流场示意图，用以反映植物随着底坡高程变化对流场的影响。

在木本植被沿斜坡均匀覆盖的情况下，波浪传播至在水平中段区域波高仍有明显的上升趋势，但随着植被拖曳力的持续作用，波浪传播至x=16.0 m后有明显的衰减趋势，波浪虽也产生波长变短，波前峰变陡的水平不对称现象，但波浪趋于平缓的向岸线方向爬升而不发生破碎。对应的，同样以左边界位置的波面高度HL=0.03 m为例生成的流场图中，相对于无植被作用情况，随着波前峰近岸趋于变陡，波峰处的水平速度和竖向速度变化不剧烈，流场状态趋于平稳。

有无植被作用情况下，波浪沿斜坡向上传播的不同位置处的波峰下沿垂线上的水平速度矢量分布如图12所示。波浪在沿斜坡爬升过程中，虽在前半段水平速度均随着爬坡而增大，但在木本植被的消浪作用下，波浪的波面高度趋于稳定，波峰下沿垂线上的水平速度值相对较小，尤其在中后段a、b两种情况的波峰垂线上的水平速度分布差异明显：随着植被消能作用的叠加，有植被作用下的波浪水平速度持续减小，趋于平缓的向坡顶爬升；而无植被消能作用下的波浪在近坡顶处的水平速度仍持续增大，波浪波动较大，趋于破碎。

图8 濒临破碎波峰下水平速度垂向分布比较Fig.8 Comparison of vertical distribution of horizontal velocity under breaking wave crest field

图9 二阶Stokes波模拟结果验证Fig.9 Verification of second-order Stokes wave simulation results

图10 植被作用下波浪沿斜坡（i=1/30）传播波形Fig.10 The wave form diagram of the wave propagation along slope (i=1/30) with vegetation

图11 植被作用下波浪沿斜坡（i=1/30）传播流场Fig.11 The flow field diagram of the wave propagation along slope (i=1/30) with vegetation

图12 无植被作用（a）和有植被作用（b）波峰下水平速度垂向分布对比Fig.12 Comparison of vertical distribution of horizontal velocity under wave crest field with vegetation (a) and without vegetation (b)

图13 无植被作用（a）和有植被作用（b）波峰位置波面高度沿程变化Fig.13 The variation of wave surface height on wave crest along distance with vegetation (a) and without vegetation (b)

为更直观的体现木本植被的消浪作用，给出了有无植被两种情况下波浪的波面高度沿程变化（图13），且考虑波浪爬升的水面变化和二阶stokes波的波形特点，以沿程波峰位置处的波面高度代替波高进行分析。由图可见，植被作用下的波面位置变化的总体趋势为前半段稳定增高，而后半段稳定减小，且趋于平缓的爬升而不发生破碎。

4.2 植被特性对消浪效果的影响规律分析

考虑植被作用情况下，为综合考虑木本植被对波能衰减的影响机理，需考虑的木本植被特性为：植被密度（单位长度内植株的数量）、树枝密度（每棵植株的树枝数）、树干高度，树干直径、树枝直径和树枝倾斜度等因素。其中，由于植物在栽种后的生长过程中其树干直径和树枝直径短期内变化缓慢，因此本文不对其进行讨论，而通过分析沿程波峰位置处的波面位置变化，着重讨论植被密度、树枝密度、树干高度和树枝倾斜度对波能衰减的影响机理。本节给定的波浪参数同上节，表1为考虑不同的植被特性条件下的数值模拟工况。

4.2.1 植被密度对消浪效果的影响

图14是对应工况1～4的沿程波峰位置处的波面位置变化，x轴为波浪向岸传播的水平方向距离，z轴为波峰位置处的波面位置，图中对应各个工况折线沿程的最后数据点为波浪传播至破碎前的波峰位置处的波面高度值。

表1 植被特性模拟工况表Table 1 The arrangement of simulation condition about vegetation characteristics

图14 植被密度与波面位置关系Fig.14 The relationship between vegetation density and wave surface position

植被作用下，工况1～4通过控制树干高度、树枝倾斜度为定值，即通过同比例改变树干、树枝的数目，实现植被分布密度的改变。树干与树枝的比例为1:6，单位长度内的植株的变化范围为 0 .1≤NL≤3，其中NL=0.1（工况1）的情况为接近无植被作用的极端情况。在这一范围内，工况1与工况2、3、4的沿程波面高度变化对比明显，且沿程高度差不断增大。在x∈[22, 24]这一区段内，波面高度工况1与工况4的高差值达到了0.13 m。而工况2～4的波面高度均先呈现稳定的、较小幅度的爬升，继而在水平传播方向上的中后段出现较为明显的下降趋势。产生这种现象可能主要由于以下两个原因：一方面，在植被消浪和倾斜底坡两种条件综合影响的情况下，波浪存在沿程消能波高减小和斜坡爬升波高增大两种趋势，在计算区域的中前段，倾斜底坡对波浪爬升的影响作用较明显，而在中后段，植被的沿程拖曳力对波能衰减的作用逐渐显现，并出现较大幅度的波面高度的减小；另一方面，当水位较低时，即波浪传播至近岸水深较浅的区域时，此时波能较小，且树干对波浪施加与其传播方向一致的拖曳力，起主要消浪作用，使波能迅速衰减。

对比工况2、3、4的情况，工况3的植被密度为工况2的2倍，工况4的植被密度为工况2的3倍。受植被密度影响，工况3的沿程波面高度略小于工况2，且波浪破碎区域相同；工况4虽然植被密度较大，但在沿程x＜26 m的区段内，工况4相对于工况2、3波面高度较高，而在后段波面高度减弱更加明显，波浪更加平稳地向近岸方向传播而不发生破碎。工况4这种沿程波面高度的变化可能是由于，在植被覆盖密度较大的情况下，植物对水体的拖曳力叠加，植被区域内水体波动较为剧烈，导致了在x＜26 m的区段内波面位置较高的现象，而与此同时，也增加了能量的损耗，使透过植被区后的波能大大减弱，对应的，波面位置则出现明显的减弱现象。由此可见，单纯增加植被的覆盖密度来增加其对波浪的衰减作用的方式并不合理，一方面经济性欠佳，另一方面此方式也可能会增加沿程水体的波动，而对生态、岸滩防护等起到不利的影响。因此，在本小节的4种工况中，工况2更符合自然植被的实际分布情况，也较为经济可行，综合效果更好。

4.2.2 树枝密度对消浪效果的影响

相对于工况2而言，工况5～7对树枝密度进行了调整。图15是对应工况2、5～7的沿程波峰位置处的波面位置变化，图中坐标轴和折线的表示意义与图14一致。由图15可以看出，在x∈[0, 12]这一区段内，随着树枝密度的增大，沿程波面高度的变化不明显，而后则出现波面高度的较大差异：树枝密度越大，波面位置的变化越趋于平缓。在x∈[0, 12]这一区段内，植被覆盖段区域为10 m，其长度约等于区域内的波长，而植被对波浪消能的作用是在传播中不断叠加的，而在沿程第一段波长的长度范围内，树枝的效果尚未体现。但随着波浪不断向前传播，树枝密度对波面位置的影响愈发显著，其中工况2的沿程波面高度最为稳定，且在x∈[18, 20]这一区段内，工况2与其他工况相比有较明显的波浪衰减，波面高差约为0.05 m。

图15 树枝密度与波面位置关系Fig.15 The relationship between branch density and wave surface position

4.2.3 树干高度对消浪效果的影响

工况8～11对树干高度进行调整，即改变了水体的上、下分层界限的高度（对应图3中的d1）。图16是对应工况2、8～11的沿程波峰位置处的波面位置变化，图中坐标轴和折线的表示意义与图14一致。类似的，在沿程第一段波长的区段内，沿程波面高度的区别不大，树干高度的变化对波浪衰减的影响体现不明显，而在此之后的区段，树干高度对波面高度变化的影响效果较为明显。

图16 树干高度与波面位置关系Fig.16 The relationship between vegetation density and wave surface position

树干高度的变化改变了图3中水体上下层分界线的位置，依据波浪传播过程中大部分波浪能量集中在表层水体内这一特点，本文讨论树干高度变化的影响主要分析其对表层水体的影响，其主要体现在以下两个方面：（1）改变树干和树枝在水体中的作用区域，即改变了水体不同区域的拖曳力大小；（2）树干-树枝组合作用、树枝主导作用、树干主导作用3种对水体表面紊动消能的影响。

由图16可知，植被的树干高度过高或过低，其工况对应的沿程波浪的衰减效果均较差：工况8和工况11分别对应树干高度hL=0.2 m和hL=1.0 m的两种情况，相对其他3种工况（工况2、9、10），其沿程波面高度值总体较大。其中，工况8对应的波浪在沿程传播过程中，在区域I和区域II的中前段，即水深较深的区域，表面波完全受树枝的拖曳力作用。而由于树枝的直径较小、作用力方向与波浪传播方向有倾角等原因，其波浪衰减效果较差。工况11对应的波浪在区域I同时受到树干和树枝的拖曳力作用，随着波浪沿斜坡向上爬升，水深逐渐减小，树干也逐渐完全代替树枝对波浪施加与其传播方向一致的拖曳力。

对比5种工况，分析其波浪沿程波面高度的变化规律可以得到结论：在水深较深、波能较大的情况下，相比于树枝主导作用和树干主导作用两种情况，树枝与树干组合作用于水体表面时，消浪效果更好。这可能是由于组合作用下，既有树干较强的拖曳力作用，也有树枝形态对水体的扰动作用，增加了水体紊动造成的波能损耗。而对比工况8、9与其他3种工况在中后段的波面高度变化情况，可以发现：在水深较浅、波能衰减程度较大的情况下，若树干代替树枝对波浪施加拖曳力，波能衰减较明显，树干起主要消浪作用。综合以上几种工况，工况9（树干高度hL=0.4 m）消浪效果最好。

4.2.4 树枝倾斜角度对消浪效果的影响

调整树枝的倾斜角度，即改变树枝与波面线的夹角。当其夹角为90°时，即树枝对波浪施加与其速度方向一致的拖曳力，消浪效果最好。由图17，对比4种工况，可以发现树枝平均角度 θ =π/6（工况9）和θ=π/3（工况14）两种情况的消浪效果较优。在波浪传播过程中，波形存在峰谷不对称和波峰的前后不对称，且随着波浪不断向近岸方向传播，波形的不对称现象也在逐渐加剧。因此可以认为工况9和工况14在波浪的传播过程中，相对于其他两种工况，树枝作用于波浪的拖曳力方向与波浪速度方向总体更趋于一致，能增加入射波浪在竖直方向上势能和水平方向动能的消耗，消浪效果较优。

图17 树枝角度与波面位置关系Fig.17 The relationship between branch tilt angle and wave surface position

4.3 波浪因素对消浪效果的影响规律分析

采用无因次参数分析影响植被消浪效果的波浪因素，波浪参数一般包括：相对波高Hs/d，波陡H/Ls（ δ），相对水深d/Ls（表2），植被覆盖区域取消浪效果相对较优的工况9的植被特性参数。其中，植被覆盖区域植被特性不变时，相对水深改变的同时会引起树干、树枝入水深度的改变，并产生影响，其本质上与树干高度对消浪效果的影响机理相同，故本节主要讨论相对波高Hs/d和波陡H/Ls（δ）两个波浪因素对消浪效果的影响规律。

由于波浪因素已发生改变，故单纯分析沿程波面位置的变化不具有代表性，故取有植被作用下的沿程波高为HT，无植被作用下沿程波高为HI，采用透浪系数Kt=HT/HI研究植被作用下对波高的影响，透浪系数越小，消浪效果越好。由3.1节的分析可知，随着水深变浅，波长沿程逐步减小，故将沿程分为[0, 7)、[7, 14)、[14, 21)、[21, 27)、[27, 32)5 个不等间距的区间段，并对区段内的透浪系数及其对应的沿程波高值的变化规律进行分析。

表2 植被特性模拟工况表Table 2 The arrangement of simulation condition about vegetation characteristics

4.3.1 相对波高HS/d对消浪效果的影响

图18是对应工况a～d的透浪系数随相对波高变化，x轴为波浪向岸传播的水平方向距离，z轴为透浪系数大小。图19是相应的沿程波高随相对波高变化图，z轴为波高大小。图中折线上最后一个数据点均为波浪传播至破碎前的波峰位置处的波面高度值。相对波高的变化范围为 0 .10≤H0/d≤0.25，透浪系数沿程变化规律可分为两段：在区域I和区域II的前段，波浪爬升、水流的惯性等因素导致波浪的变形情况较复杂，透浪系数无明显规律；而在区域II的中后段，波浪变形趋于稳定，在这一范围内，随着相对波高的增大，透浪系数减小趋势明显。这是由于相对波高越大，枝干对水体表层作用的范围越大，波能损耗越大。对应的，由图19可见，在区域I和区域II的前段，沿程波高值大小与相对波高大小呈正比关系，而随着沿程趋于近岸方向，波高值的差距呈缩小趋势，虽然波高值的变化还受到爬升、破碎的影响，但总体而言，随着相对波高的增大，沿程波高值衰减速度相对较快。

图18 相对波高与透浪系数关系Fig.18 The relationship between relatively wave height and transmitted coefficient

图19 相对波高与沿程波高关系Fig.19 The relationship between relatively wave height and wave height along the path

4.3.2 波陡HS/LS(δ)对消浪效果的影响

图20是对应工况a、e～g的透浪系数随波陡变化，x轴为波浪向岸传播的水平方向距离，z轴为透浪系数。图21是相应的沿程波高随波陡变化，z轴为波高大小。

图20 波陡与透浪系数关系Fig.20 The relationship between wave steepness and transmitted coefficient

图21 波陡与沿程波高关系Fig.21 The relationship between wave steepness and wave height along the path

波陡的变化范围为 0 .029≤δ≤0.048，透浪系数沿程变化规律可分为两段：在区域I和区域II的前段，与图18类似的，波浪爬升、水流的惯性等因素导致波浪变形变化较复杂，透浪系数无明显规律；而在区域II的中后段，波浪变形趋于稳定，在这一范围内随着波陡的增大透浪系数呈现增大趋势。由图21可以看出，在给定的波陡范围内，植被的消浪效果较好，波浪沿程平缓的向岸方向传播，波浪沿程不会发生破碎等剧烈紊动现象，即随着波陡的增大，不会产生的较大的波能损耗。而此时，波长对消浪效果起主导作用的可能性较大，大波浪较大的往复运动幅度增加了水体与植被之间的相互作用，即大波浪作用时的消浪效果要好于小波浪的消浪效果。类似的，杨建民[18]、郭晓宇等[26]根据研究植被消浪的试验数据及现象观察也得到了相同的结论。根据线性波理论，同一波高条件下，长波的水质点速度沿水深变化较短波慢，相应的长波受到的影响可能会比较显著[26]。

5 结论

针对近岸种植木本植物这一新型的海岸生态软防护措施，本文考虑植被枝干的分区作用，在N-S方程中分别考虑树枝和树干的拖曳力影响，提出了木本植被保护下波浪沿斜坡爬升的表面波衰减理论模型，并采用有限差分法离散方程和MAC方法来跟踪自由曲面上的水颗粒轨迹。以波浪沿1/30的斜坡爬升为算例，对比讨论了有无植被作用下波浪沿倾斜海滩上传播过程，并将算例结果与以往试验结果规律进行对照，验证了数值模型的有效性。在木本植被均匀覆盖岸坡这一条件下，分别讨论了植被分布密度、单株植物树枝密度、树干高度、树枝倾斜角度4种植被特性和相对波高、波陡两种波浪因素对植被消浪效果的影响，得到植被消浪的基本规律。

针对植被特性对消浪效果的影响规律研究，结果表明：（1）单位长度内的植株的变化范围为0.1≤NL≤3时，在保证消浪效果的情况下，从经济性和岸滩防护的两方面考虑，NL=1时效果最好，实际工程中，单纯增加植被的覆盖密度来增加其对波浪的衰减作用的做法有待进一步讨论，其可能导致植被区域内水体波动较为剧烈，局部区段内波面高度变化较大，不利于生态和岸滩防护；（2）单植株的树枝数变化范围为3≤NU≤6时，树枝密度越大，波面位置变化越趋于平缓，在NU=6时消浪效果最好；（3）树干高度变化范围为0.2 m≤hL≤1.0 m时，树干高度为hL=0.4 m消浪效果最好，在水深较深，波能较大的情况下，树枝-树干组合情况增加了水体紊动造成的波能损耗，消浪效果较好；在水深较浅、波能衰减程度较大的情况下，树干起主要消浪作用；（4）树枝倾斜角度在 π/6≤ θ≤π/3范围内时，θ=π/6和 θ=π/3两种情况树枝作用于波浪的拖曳力方向与波浪速度方向总体更趋于一致，消浪效果较优。针对波浪因素对消浪效果的影响规律研究，结果表明：相对波高的变化范围为 0.10≤HS/d≤0.25，透射系数总体随着波陡与相对波高的增加呈缓慢减小趋势；波陡变化范围为 0 .012≤δ≤0.048，在这一范围内，波浪沿程传播平缓，波陡的变化对消浪效果影响不明显，波长对消浪效果起主导作用，大波浪作用时的消浪效果好于小波浪。文中的计算结果和消浪的基本规律可为实际的护岸工程和生态景观设计提供参考。

植被护岸这一软防护措施在实际使用中，多设置在倾斜底坡的海滩或近岸结构上，而由于波浪与植被本身作用十分复杂，因此，提出考虑植被形态和斜坡爬升的表面波衰减理论模型、分析波浪作用下植被参数与波要素对消浪效果的影响规律，对工程实际具有重要的意义。但由于波浪在植被覆盖区域围内的紊动情况复杂，且未考虑树枝部分的柔性结构对消浪效果的影响，对植被特性的讨论还有待进一步研究与分析。