植被对波浪作用下床面切应力影响的数值模拟分析

葛昭佩，唐军*，赵楚嫣

( 1. 大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连 116023)

1 引言

波浪作用下的床面切应力对海床泥沙起动具有关键性作用，是模拟分析近岸泥沙输移、岸滩演变的重要参数。近岸植被具有防浪固滩护岸的作用，研究近岸植被对波浪作用下床面切应力的影响有助于理解其防护海岸的作用机制。

目前关于波浪作用下床面切应力的研究取得了一定成果。Jonsson[1]从理论公式出发，推导了层流边界层的床面切应力及边界层厚度公式，并引入波浪底摩阻系数及边界层流态判别方法，给出了紊流边界层的波浪摩阻系数经验公式；孔令双等[2]基于边界层理论，导出了波流环境中床面切应力，并建立了泥沙起动及水体挟沙力公式；蔡翠苏[3]基于水槽试验，探讨了规则波和不规则波作用下底摩阻系数计算方法；齐富康等[4]基于渤海海峡观测数据，分别应用Soulsby和Van Rijn 的模型计算波流环境中床面切应力；Lin和Zhang[5]建立三维数值模型，模拟计算线性波、Stokes 波、孤立波等条件下层流边界层床面切应力，并分别与解析解比较分析；滕涌等[6]基于ECOMSED波流耦合底边界层模型，分析了波高、近底流速及水深对波流条件下床面切应力的影响；Larsen 和Fuhrman[7–8]基于waves2Foam 模拟海啸波的传播及爬升过程，着重讨论了海啸波作用下的边界层及床面切应力并提出了一种预测海啸波作用下瞬时边界层厚度和床面切应力的方法。相对裸床，植被影响下波浪水动力特性复杂，目前关于植被水域波浪作用下床面切应力的计算以经验公式为主。Wang 等[9]基于水槽试验分析了水流和波流条件下植被水域平均流速及湍流动能，并由湍流动能推算植被水域床面切应力，发现纯水流时植被水域床面切应力沿程降低，而波流条件下植被水域床面切应力增大并在离开植被区后显著大幅降低；Reidenbach 和Thomas[10]观测了弗吉尼亚海岸保护区水动力，并应用经验公式分别计算波、流的床面切应力，发现与纯水流时相比波流条件下床面切应力大幅增大，位于大叶藻水域测点处的床面切应力始终小于裸床处且小于泥沙起动临界切应力，有效保护了海床免受侵蚀；陈家贵和沈小雄[11]基于Jonsson[1]提出的层流边界层公式，结合水槽试验数据，分析了入射波高及柔性植被对最大床面切应力的影响；李勰等[12]基于Jonsson[1]提出的边界层最大切应力公式及Luhar 等[13]提出的植被带中边界层流速公式，分析了刚性植被根茎对床面切应力沿程变化的影响。

总体来看，目前关于植被对波浪或波流作用下床面切应力影响的研究主要基于水槽试验或现场观测数据并结合经验公式分析计算，但由于植被影响下的流速剖面及雷诺应力分布复杂，经验公式未能很好地描述波浪条件下植被水域床面切应力[9,14–15]。本文以OpenFOAM 中的waves2Foam 求解器[16]为基础，在其中引入植被源项，建立植被水域波流数值水槽，通过计算全水深速度场及边界层内速度梯度而直接给出床面切应力分布。在应用已有实验数据及理论公式对数值模型验证的基础上，着重分析了入射波高、植被密度、植被淹没高度及水流对植被水域波浪作用下床面切应力的影响特性。

2 数值模型

waves2Foam 求解器是由Jacobsen 等[16]基于Open-FOAM 开发的一个用于模拟波浪的第三方模块库，提供了水流、线性波、Stokes 波、孤立波等造波文件，亦可自定义造波文件实现其他类型波浪模拟。

2.1 控制方程

waves2Foam 求解器控制方程基于OpenFOAM 标准，采用气液两相流方程。为考虑波浪与植被相互作用，在动量方程和湍流模型中引入植被源项，即拖曳阻力项和惯性力项[17]（称之为宏观数值模型，以下未特别说明的均指该数值模型）。控制方程为

传统湍流模型应用于气液两相流时，由于气液界面处较大的密度梯度[19]及湍流动能的过度评估[20]会造成大波陡波传播过程中的波高非物理衰减。为抑制此种现象，本文采用Larsen 和Fuhrman[20]提出的修正k-ε湍流模型。

表1 经验系数取值Table 1 Default values for the closure coefficient

waves2Foam 为两相流求解器，其定义体积分数α （ 0 ≤α ≤1）来描述每个网格中各相的体积比重。各相的运动满足的控制方程为

式（8）第三项是用于保持气液界面清晰的人工压缩项；ur,i为 压缩速度；压缩系数Cα∈[0，1]，其默认值为1。

2.2 植被拖曳力系数

植被拖曳力系数CD是对植被水域水动力特性进行宏观数值模拟时的关键参数，与波浪、水流、植被密切相关，直接影响数值模拟的准确性，但其取值尚无普适方法。本文将上述控制方程中的植被源项去除（令CD、CM、Ckp、Cεp4 项经验系数均取为0），通过直接考虑刚性植被外形对波流的扰动，详细模拟纯波及波流条件下植被水域的水动力特征（称之为精细数值模型，在本文中仅用于计算各工况代表CD值），通过提取植被上的波流力及其所在截面的沿植被高度平均流速并使用Morison 公式（式（10））直接计算植被拖曳力系数。

式中，FD为 植株上的拖曳力；FI为植株上的惯性力；F为植株上的波流力，Dalrymple 等[22]指出阻力与流速之间无相位差，阻力的做功可由阻力和流速的乘积在周期内积分获得，然而惯性力与流速之间存在 π/2的相位差，使得惯性力在周期内的做功为0，故上式计算过程中可忽略惯性力[18]，F根据植被上压力及附近流场由式（11）计算；u为植株所在截面的平均流速；ρ为水的密度；hv为 植被淹没高度；bv为植被直径。

2.3 边界条件

数值水槽入口与出口采用松弛区造波与消波，入口处速度场、压力场和自由波面由波浪理论公式给出，出口采用自由出流条件；床面采用无滑移边界条件；两侧壁采用周期性边界条件，该边界条件能够消除两侧壁面的影响，可将水槽简化为较窄的数值水槽以提高计算效率，且已被成功应用于波浪的数值模拟[23–24]；顶端采用压力出口边界。时间步长由库朗数自动调节，为保证计算稳定，最大库朗数设置为0.25。

3 模型验证

由于缺乏波浪作用下植被水域床面切应力实测数据，本文将分别验证植被水域波面演化和无植被时波浪作用下床面切应力。采用海岸和近海工程国家重点实验室开展的植被影响下波浪传播物理水槽试验结果[25]验证本文模型模拟植被水域波浪传播的有效性；采用理论公式验证层流边界层时的床面切应力，采用徐华等[26]的试验验证紊流边界层时的床面切应力。

3.1 植被水域波浪传播验证

选取表2 中的两种试验工况验证模型模拟植被水域波浪传播的有效性，其中CD值由精细数值模型计算，CD值计算过程以表2 中工况1 为例，数值模拟中设置数值水槽长为18 m，高为0.6 m，宽为0.12 m，取植被水域长为1.08 m，其中以圆柱代表刚性植被，整体网格尺寸 Δx=Δy=0 .03 m， Δz=0.01 m，为更好地贴合植被外形，分别加密近植被域及单株植被附近网格。边界条件同2.3 节，此外，植被表面设为无滑移边界条件。

表2 植被及波浪参数Table 2 Parameters of vegetation and waves

图1 为精细模拟中工况1 的波面演化和x=0.6 m（x=0 为植被域起点）处植被上的波浪力及截面平均速度，使用式（10）计算该处周期平均CD值，并利用同样的方法计算每一排植株处的周期平均CD值后，做平均得出空间–周期平均CD值作为该工况的代表CD值。

图1 考虑植被外形扰动的模拟结果Fig. 1 Simulation results considering the disturbance of vegetation shape

将利用精细模型计算的各工况代表CD值代入宏观模型中模拟计算。波浪经过植被水域的波面演化及流速衰减的模拟值与试验值对比如图2 和图3 所示，由图可知，模型计算结果和试验数据[25]符合性较好，表明利用精细数值模型计算CD值准确并且使用宏观数值模型能够准确模拟植被水域波浪及波流水动力变化。

图2 植被水域波面演化Fig. 2 Free surface evolution along the vegetation zones

图3 植被水域流速衰减验证Fig. 3 Verification of velocity attenuation in vegetation zones

3.2 波浪作用下床面切应力验证

本文数值模型通过求解雷诺平均Navier-Stokes 方程得出波浪作用下全水深速度场及边界层内速度梯度，并由公式（12）计算床面切应力。

为验证模型计算床面切应力的准确性，此处设置3 组验证工况，其中工况1、工况2 为本试验工况，工况3 为徐华等[26]的试验工况，波浪参数如表3 所示。工况1 和工况2 的边界层雷诺数均小于1.26×104，依据Jonsson[1]提出的判定标准，两者均属于层流边界层，工况3 则属于光滑紊流边界层。

表3 验证工况参数Table 3 Parameters of verification conditions

边界层流速剖面及床面切应力验证结果如图4至图6 所示。其中，工况1 采用线性波理论解验证，工况2 采用五阶Stokes 波理论解验证，工况3 采用试验测量值[26]验证。从图中可以看出，待波浪场稳定后，工况1 模拟值与理论值符合性较好；工况2 床面切应力在波谷时模拟值较理论值偏大，但差值小于10%；工况3 模拟值与试验测量值在一个完整周期内总体符合良好，表明该数值模型能够准确计算波浪作用下床面切应力。

图4 理论值与模拟值对比（工况1）Fig. 4 Comparison between theoretical and simulated values (case 1)

图5 理论值与模拟值对比（工况2）Fig. 5 Comparison between theoretical and simulated values (case 2)

图6 床面切应力对比（工况3）Fig. 6 Comparison of bed shear stress (case 3)

4 模拟结果分析

为分析不同入射波高、植被密度、植被淹没高度及水流流速条件下，植被水域床面切应力的分布特征，本文共模拟计算22 种工况，见表4。其中工况1、2、3、12、14、15、18、19、20 为文献[25]的试验工况，这9 种工况的波面演化及流速衰减模拟值与试验值均符合良好。

表4 模拟工况参数Table 4 Parameters of numerical simulation conditions

纯波时，波浪经过植被水域后波高衰减的同时床面切应力也随之减小，见图7（x=0 m 处为植被起点），这是由于植被的阻水作用，近底流速沿植被水域衰减使得边界层流速梯度减小所致（图8）。图7 中植被前后的波面与床面切应力均存在 π/4的相位差，说明植被对波高和床面切应力的衰减作用是同步的。此外，该相位差与式（13）、式（14）相一致，进一步验证了该数值模型的准确性。

图7 波面及床面切应力（工况1）Fig. 7 Free surface and bed shear stress (case 1)

图8 近底流速剖面（工况1）Fig. 8 The near-bottom velocity profile (case 1)

当植被处于淹没状态时，植被的阻水效应仅发生在冠层以下，使这一区域流速减小且阻水作用随着淹没高度的增大而增大，而冠层以上由于过水断面的减小使得流速增大且超过无植被时的流速，故在植被顶端形成较大的速度梯度（图9）；此外，由于边界层的作用在近底处存在一速度峰值。上述水平速度沿水深的分布规律与Liu 等[29]的试验结果一致。

图9 淹没植被沿水深速度剖面Fig. 9 Longitudinal velocity profile of submerged vegetation

同向水流的存在会增加波浪水质点速度的非线性，同时增大正向床面切应力幅值，减小负向床面切应力幅值，使得一个波周期内正向切应力时间延长，负向切应力时间缩短，如图10 所示。

图10 不同流速下床面切应力变化Fig. 10 Variation of bed shear stress under different current velocity

5 结论

本文OpenFOAM 中的waves2Foam 求解器建立了含植被水域的波流数值水槽，模拟分析了波浪、水流、植被要素对植被水域床面切应力衰减的影响。模拟结果表明：当波浪传播到植被前端时，床面切应力会出现增大现象，并随着入射波高增大、水流流速增大、植被密度增大此种现象愈加明显；纯波时，波高和植被密度的增大均会导致植被水域床面切应力衰减幅度的增大，并且床面切应力在植被水域前段的衰减幅度较大，后段衰减幅度减小；对于完全淹没植向床面切应力幅值减小，使得一个波周期内正向切应力时间延长负向切应力时间缩短；当水流流速较小时，床面切应力的衰减率与无水流时相近，之后随着流速增大衰减率随之增大。

图11 不同工况下植被水域最大床面切应力分布Fig. 11 Distribution of maximum bed shear stress in vegetation zones under different conditions

图12 不同流速下最大床面切应力衰减率Fig. 12 Decay rate of maximum bed shear stress at different current velocities