刚性沉水植物影响下波浪传播及泥沙悬浮数值模拟

摘要

由于植物引起的紊动作用，含刚性植物床面中悬浮泥沙浓度与无植物裸床相比显著增加，而通过平均流速计算底床切应力的传统泥沙模型无法模拟出该现象。因此，在含植物水流水沙运动物理模型试验的基础上，构建了基于Flow–3D的含刚性沉水植物条件下波浪传播的三维数学模型，模拟了植物影响下波浪动力特征和泥沙悬浮过程，同时从紊动能角度修正希尔兹数，对泥沙模块进行了改进。与实测数据相比，该模型可较精确地模拟出植物引起的整体水流流速减小和局部冠层顶部处的流速增大、水体紊动增强以及紊动能在波周期内出现两个峰值的现象。与原始泥沙模块相比，改进的模型考虑了植物尾流紊动对泥沙运动的影响，可提高含植物床面泥沙悬浮模拟的精度。

关键词

Flow–3D数值模拟; 刚性沉水植物; 波浪流速; 紊动能; 泥沙模块改进; 泥沙悬浮

海岸带、滨海湿地等区域普遍存在的水生植物，对水动力和水环境具有重要影响。一方面，植物可以减小水流流速，削弱波浪对岸滩的冲击；另一方面，植物也会改变水体的紊动特征，影响泥沙输运和污染物迁移过程。因此，研究含植物水流水动力特征及泥沙输运过程对岸滩保护和水环境修复具有重要意义。

目前，大多数含植物水流数学模型将植物群设置为多体结构或多孔介质进行模拟。多体结构模型是将每根植物视为独立的个体，进而可将整体植物群等效为柱体群结构，利用网格读取每根植物的边界。如Chang等^［1］利用分离涡模拟（DES），构建了圆柱群阵列，对每根植物的几何形状进行了精细模拟，分析了淹没刚性植物对水流流速和紊动的影响，在植物尾迹区观察到U形涡旋。多孔介质模型是用一个带有孔隙率的整体代替植物群，通过在计算方程中添加概化固体阻力源项来模拟。如Zhan等^［2］采用多孔介质模拟方法，通过设置统一的孔隙率将植物群概化为等间距圆柱群，并引入一个非恒定的惯性阻力系数来模拟不均匀的植物阻力分布。与多体结构模型相比，多孔介质模型虽然可以节省大量计算时间，但可能存在植物群内部水体结构模拟失真的问题^［3］。如Yu等^［4］利用Fluent模型对比了两种方法下刚性挺水植物与水流相互作用，发现多体结构模型可以相对较好地模拟出植物影响下的平均流速场和紊流场，多孔介质模型只能模拟植物群整体产生的大尺寸紊动而不能很好地模拟单根植物引起的的小尺寸紊动。因此，为了更加精确地模拟植物对水动力特征的影响，获得植物群内部的水流结构和泥沙运动规律，本研究选用多体结构方法模拟刚性沉水植物。

近年来越来越多的研究发现，植物的存在使得水体中悬浮泥沙浓度显著增加，含植物床面中植物引起的紊动对泥沙起动及悬浮起主导作用^［5-6］。因此，仅考虑平均流速的床面切应力并不适用于含植物水流中泥沙起动与输运的模拟。Tinoco和Coco^［5， 7］用紊动强度代替平均流速，采用物理模型试验数据，对传统临界希尔兹（Shields）数进行修正，提出了适用于波浪、水流条件下的含植物床面泥沙起动判断方法，但该方法是基于特定的物理模型试验结果进行的推导，其应用范围仍需进一步探讨。Lou等^［8］进一步修正Shields数并应用于NHWAVE模型，模拟了植物影响下的泥沙悬浮过程，并获得了较为精确的结果。然而，上述改进的模型需要对较多参数进行校正，其中一些参数的物理意义难以得到合理解释。因此，本文在前期物理模型试验的基础上，基于当前被广泛认可的利用流速和植物特征计算植物尾流紊动的方法^［9-10］，从理论上采用紊动能对传统Shields数进行修正，从而获得可应用于含植物床面的泥沙起动与悬浮的预测方法。该方法具有较强的物理意义，且并未引入新的参数。应用该方法，改进了Flow–3D模型中的泥沙模块，开展了波浪通过刚性沉水植物区时泥沙悬浮的模拟。

1 数学模型建立

1.1　基本方程

1.1.1　水流方程

Flow–3D模型中笛卡尔坐标系下的控制方程如下：

连续性方程：

$$V_{\mathrm{F}}\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_i}(\rho u_i{A}_i)=\mathrm{0}$$

（1）

动量方程：

$$\frac{\partial u_i}{\partial t}+\frac{\mathrm{1}}{V_{\mathrm{F}}}{u}_j{A}_j\frac{\partial u_i}{\partial x_j}=-\frac{\mathrm{1}}{\rho }\frac{\partial p}{\partial x_i}+G_i+f_i$$

（2）

式（1）—（2）中：V_F为流体的体积分数；ρ为流体密度；t为时间；x_i为笛卡尔坐标（x，y，z）；u_i为速度分量；A_i为i向流体的面积分数；p为压力；G_i为i向体加速度；f_i为i向粘滞力加速度。

采用Yakhot和Orszag^［11］提出的RNG k-ε模型作为紊流模型。紊动能k_T和紊动耗散率ε_T满足如下方程：

$$\frac{\partial k_{\mathrm{T}}}{\partial t}+\frac{\mathrm{1}}{V_{\mathrm{F}}}{u}_i{A}_i\frac{\partial k_{\mathrm{T}}}{\partial x_i}=P_{\mathrm{T}}+G_{\mathrm{T}}+D_k-{\epsilon }_{\mathrm{T}}$$

（3）

$$\frac{\partial {\epsilon }_{\mathrm{T}}}{\partial t}+\frac{\mathrm{1}}{V_{\mathrm{F}}}{u}_i{A}_i\frac{\partial {\epsilon }_{\mathrm{T}}}{\partial x_i}=\frac{{C_{\mathrm{1}}}_{\epsilon }\cdot {\epsilon }_{\mathrm{T}}}{k_{\mathrm{T}}}(P_{\mathrm{T}}+{C_{\mathrm{3}}}_{\epsilon }\cdot G_{\mathrm{T}})+D_{\epsilon }-{C_{\mathrm{2}}}_{\epsilon }\frac{{\epsilon }_{\mathrm{T}}^{\mathrm{2}}}{k_{\mathrm{T}}}$$

（4）

式（3）—（4）中：P_T为速度梯度引起的紊动能产生项；G_T为浮力引起的紊动能产生项；D_k和D_ε为扩散项；C_1ε、C_2ε和C_3ε为模型参数。

1.1.2　泥沙运动方程

Flow–3D模型中，以临界切应力作为判断泥沙起动的条件，采用Soulsby-Whitehouse公式^［12］计算量纲一化临界Shields数θ_cr：

$${\theta }_{\mathrm{c}\mathrm{r}}=\frac{\mathrm{0.3}}{\mathrm{1}+\mathrm{1.2}{d}_{\mathrm{*}}}+\mathrm{0.055}\left[\mathrm{1}-\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{0.02}{d}_{\mathrm{*}})\right]$$

（5）

$$d_{\mathrm{*}}=d_{\mathrm{s}}{\left[\frac{\rho ({\rho }_{\mathrm{s}}-\rho )‖g‖}{{\mu }^{\mathrm{2}}}\right]}^{\mathrm{1}/\mathrm{3}}$$

（6）

式中：d_*为量纲一化泥沙粒径；d_s为泥沙粒径；ρ_s为泥沙密度；||g||为重力加速度的大小；μ为流体动力粘度。

泥沙挟带上升速度u_lift用于计算床沙转化为悬浮状态的量，采用Mastbergen和Berg^［13］的理论公式：

$$u_{\mathrm{l}\mathrm{i}\mathrm{f}\mathrm{t}}=\alpha n_{\mathrm{s}}{d}_{\mathrm{*}}^{\mathrm{0.3}}(\theta -{\theta }_{\mathrm{c}\mathrm{r}}{)}^{\mathrm{1.5}}\sqrt[]{\frac{‖g‖d_{\mathrm{s}}({\rho }_{\mathrm{s}}-\rho )}{\rho }}$$

（7）

$$\theta =\frac{\tau }{‖g‖d_{\mathrm{s}}({\rho }_{\mathrm{s}}-\rho )}$$

（8）

式中：α为泥沙挟带系数；n_s为床面法向向量；θ为Shields数；τ为床面切应力。

悬移质泥沙扩散方程为

$$\frac{\partial C_{\mathrm{s}}}{\partial t}+\nabla \cdot \left\{\left[\overline{u}+\frac{v}{d_{\mathrm{s}}}\left({\left(\mathrm{10}.{\mathrm{36}}^{\mathrm{2}}+\mathrm{1.049}{d}_{\mathrm{*}}^{\mathrm{3}}\right)}^{\mathrm{1}/\mathrm{2}}-\mathrm{10.36}\right)\frac{g}{‖g‖}\cdot \frac{C_{\mathrm{s}}}{{\rho }_{\mathrm{s}}}\right]C_{\mathrm{s}}\right\}=\nabla \cdot \nabla \left(D{C}_{\mathrm{s}}\right)$$

（9）

式中：C_s为悬浮泥沙的质量浓度；$\overline{u}$为水沙混合物的速度，利用u_lift和水流速度相乘的通量计算得到；v为流体运动粘度；D为扩散系数。

1.2　改进的泥沙模块

由式（7）和（8）可知，Flow–3D原始泥沙模块是基于流速计算的床面切应力来判断和计算泥沙起动及悬浮的。本研究考虑植物引起的紊动作用对泥沙输运的影响，对该模型进行了改进。

对于裸床面（无植物），当单向流通过时，由床面引起的紊动k_Tb与床面剪切τ_c之间存在线性关系^［14］：

$${\tau }_{\mathrm{c}}/(\rho k_{\mathrm{T}\mathrm{b}})=\mathrm{0.19}$$

（10）

同时，τ_c也可用平均流速U和床面摩擦系数C_f直接计算^［12］：

$${\tau }_{\mathrm{c}}=\rho C_{\mathrm{f}}{U}^{\mathrm{2}}$$

（11）

由式（10）和（11）可得：

$$k_{\mathrm{T}\mathrm{b}}=C_{\mathrm{f}}{U}^{\mathrm{2}}/\mathrm{0.19}$$

（12）

Tang等^［15］基于水槽试验，初步证明了式（12）同样适用于波浪条件，得出：

$$k_{\mathrm{T}\mathrm{b}}=C_{\mathrm{w}}{U}_{\mathrm{w}}^{\mathrm{2}}$$

（13）

式中：C_w为比例因子，取值为C_w = 0.02 ± 0.01；U_w为波浪水平振荡速度幅值。

将式（13）代入Flow–3D原始泥沙模块Shields数的计算公式（8），将其改写为基于紊动能的Shields数θ₀：

$${\theta }_{\mathrm{0}}=\frac{\frac{\mathrm{1}}{\mathrm{2}}\rho f_{\mathrm{w}}{U}_{\mathrm{w}}^{\mathrm{2}}}{‖g‖d_{\mathrm{s}}({\rho }_{\mathrm{s}}-\rho )}=\frac{\frac{\mathrm{1}}{\mathrm{2}}\rho f_{\mathrm{w}}\frac{k_{\mathrm{T}}}{C_{\mathrm{w}}}}{‖g‖d_{\mathrm{s}}({\rho }_{\mathrm{s}}-\rho )}$$

（14）

式中：f_w为波浪条件下底部摩擦系数。

Tanino和Nepf^［9］根据水槽试验结果提出了单向流条件下植物尾流紊动k_Tv与水流平均流速U的关系公式：

$$k_{\mathrm{T}\mathrm{v}}={\delta }^{\mathrm{2}}{\left[C_{\mathrm{D}}\frac{l}{d}\frac{\varphi }{(\mathrm{1}-\varphi )\mathrm{\pi }/\mathrm{2}}\right]}^{\mathrm{2}/\mathrm{3}}{U}^{\mathrm{2}}$$

（15）

式中：δ为比例因子；C_D为植物拖曳力系数；ϕ为冠层体积分数，对圆柱群阵列，$\varphi =\frac{\mathrm{1}}{\mathrm{4}}n\mathrm{\pi }{d}^{\mathrm{2}}$，其中n为植物密度，d为植物直径。当n较小使得植物间距S > d/0.56时，l = d，δ = 1.1；当n较大使得S ≤ d/0.56时，l = S，δ = 0.88。Zhang等^［10］利用物理模型试验证明了该公式可应用于波浪条件，但需满足A_w/S > 1，其中A_w为波浪轨道偏移量，此时δ将减小为0.76。结合式（13）和（15），波浪条件下，含植物尾流紊动的床面总紊动能为

$$k_{\mathrm{T}}={\delta }^{\mathrm{2}}{\left[C_{\mathrm{D}}\frac{l}{d}\frac{\varphi }{(\mathrm{1}-\varphi )\mathrm{\pi }/\mathrm{2}}\right]}^{\mathrm{2}/\mathrm{3}}{U}_{\mathrm{w}}^{\mathrm{2}}+C_{\mathrm{w}}{U}_{\mathrm{w}}^{\mathrm{2}}$$

（16）

将式（16）代入式（14），可得到含植物床面的Shields数θ_ϕ为

$${\theta }_{\varphi }=\frac{\frac{\mathrm{1}}{\mathrm{2}}\rho f_{\mathrm{w}}\frac{{\delta }^{\mathrm{2}}{\left[C_{\mathrm{D}}\frac{l}{d}\frac{\varphi }{(\mathrm{1}-\varphi )\mathrm{\pi }/\mathrm{2}}\right]}^{\mathrm{2}/\mathrm{3}}{U}_{\mathrm{w}}^{\mathrm{2}}+C_{\mathrm{w}}{U}_{\mathrm{w}}^{\mathrm{2}}}{C_{\mathrm{w}}}}{‖g‖d_{\mathrm{s}}({\rho }_{\mathrm{s}}-\rho )}=\left\{\mathrm{1}+\frac{{\delta }^{\mathrm{2}}}{C_{\mathrm{w}}}{\left[C_{\mathrm{D}}\frac{l}{d}\frac{\mathrm{2}\varphi }{(\mathrm{1}-\varphi )\mathrm{\pi }}\right]}^{\mathrm{2}/\mathrm{3}}\right\}{\theta }_{\mathrm{0}}$$

（17）

将式（17）代入式（7），可得到波浪条件下含植物床面的泥沙挟带上升速度u_lift，ϕ的计算公式（18）。ϕ = 0时，θ_ϕ = θ₀。用u_lift，ϕ替换Flow–3D原始泥沙模块中泥沙挟带上升速度u_lift，即可得到改进后的泥沙模块，以用于计算含植物床面上泥沙起动及悬浮过程。

$$u_{\mathrm{l}\mathrm{i}\mathrm{f}\mathrm{t},\varphi }=\alpha n_{\mathrm{s}}{d}_{\mathrm{*}}^{\mathrm{0.3}}({\theta }_{\varphi }-{\theta }_{\mathrm{c}\mathrm{r}}{)}^{\mathrm{1.5}}\sqrt[]{\frac{‖g‖d_{\mathrm{s}}({\rho }_{\mathrm{s}}-\rho )}{\rho }}$$

（18）

1.3　模型设定

本文建立基于Flow–3D的含刚性沉水植物水槽数学模型，利用作者前期物理模型试验数据^［16-17］对数值模拟结果进行验证。计算区域布置如图1所示，长20 m（x向），宽0.2 m（y向），高1 m（z向）。水槽入口采用波浪边界（Wave），出口采用出流边界（Outflow），底面采用无滑移壁面边界（No-slip Wall），其他三面均采用对称边界（Symmetry）。

图1  计算区域布置示意图（单位：m；尺寸不按比例）

Fig. 1  Schematic of computational domain（unit: m, not to scale）

模型中采用直径为8 mm的圆柱体构建长6 m，宽0.2 m的植物区（图1），共设置3组植物布置方案（图2），其中工况a为裸床对比工况（无植物）；工况b中植物密度n为200株·m^－2，高度h_v为0.4 m；工况c中n为400株·m^－2，h_v为0.4 m；为模拟真实生态环境中植物高度的空间变化特征，在工况d中布置4种不同高度的植物（0.2、0.4、0.6、0.8 m），其密度和工况c相同。在模拟植物区底部铺设0.1 m厚的沙床，泥沙粒径为0.16 mm，密度为1.45 g·cm^－3。

图2  模拟植物分布示意图（单位：cm）

Fig. 2  Schematic of mimic vegetation canopy（unit: cm）

模型中设置4组波浪工况W1~W4，波高分别为0.06、0.08、0.10、0.12 m，波周期均为1.8 s。采用五阶斯托克斯波进行模拟，水槽末端设置海绵层用于消波。为节省计算时间、确保计算精度，本模型采用双重网格嵌套的方式，设置粗网格覆盖整个水槽区域，耦合细网格对植物区进行加密。粗网格的水平间距均为0.04 m，垂向平均间距为0.025 m，在水面和底部对网格进行局部加密。细网格垂向间距和粗网格相同，水平间距选择三种不同尺寸进行网格敏感性分析。结果表明，与间距为0.001 m网格的模拟结果相比，0.002 5 m间距网格相差小于10 %，而0.004 m间距网格相差大于20 %。兼顾计算精度和效率，选择细网格间距为0.002 5 m。

采用上述尺寸的双重网格嵌套模式，对模型的造波能力进行验证。无植物时W4条件下x = 3.5 m和x = 19.98 m处（入口边界定为x = 0）的计算值与理论值对比如图3所示。在x = 3.5 m处，可以生成稳定的波浪，平均波高约0.12 m，计算值与五阶斯托克斯波理论值及实测值均吻合较好。在x = 19.98 m处，波面升高接近于0，说明波浪在海绵层内充分衰减，不会在出口边界产生反射波影响模拟结果。因此，本模型能够准确模拟波浪传播，可用于研究植物‒波浪‒泥沙之间的相互作用。

图3  不同位置处波面升高计算值与五阶斯托克斯波理论值及实测值对比

Fig. 3  Comparison of numerical, fifth-order Stokes theoretical and measured results of wave trains at different locations

1.4　模型精度评价

本文采用均方根误差（R_MSE）和一致性指数（I_A）定量评估数值模型的模拟效果。

$$R_{\mathrm{M}\mathrm{S}\mathrm{E}}=\sqrt[]{\frac{{\displaystyle \sum _{i=\mathrm{1}}^n}(P_i-O_i{)}^{\mathrm{2}}}{n}}$$

（19）

$$I_{\mathrm{A}}=\mathrm{1}-\frac{{\displaystyle \sum _{i=\mathrm{1}}^n}(P_i-O_i{)}^{\mathrm{2}}}{{\displaystyle \sum _{i=\mathrm{1}}^n}{\left(\left|P_i-\overline{O}\right|+\left|O_i-\overline{O}\right|\right)}^{\mathrm{2}}}$$

（20）

式中：P_i为模拟值；O_i为实测值；n为样本个数；$\overline{O}$为实测值的平均值。R_MSE值越小，模拟值与实测值越接近，模拟效果越好。I_A值在0-1之间，越接近1，模拟效果越好。

2 计算结果分析

2.1　波浪动力特性分析

2.1.1　流速分析

图4为断面x = 7.5 m处波浪水平振荡速度幅值U_w（= 0.5（U_max-U_min）， U_max和U_min分别为波浪最大流速和最小流速）在各工况下的对比结果，其中h_v为工况b和c中植物高度，其值为0.4 m，水平虚线表示植物冠层顶部，后同。无植物（工况a）的4组波浪条件下，模拟值与实测值较为接近，实测与模拟的U_w均自水面向底部逐渐减小，说明该模型可以较好地模拟波浪流速（R_MSE < 0.005 m·s^－1，I_A > 0.91）。三组植物工况下（工况b~工况d），模拟值与实测值分布趋势一致，但在个别条件下出现了一定程度的误差（R_MSE < 0.014 m·s^－1，I_A > 0.75），如工况c-W1和c-W4中，冠层下方模拟值与实测值存在一定偏差，这可能是由于物理模型试验中测点位置偏差引起的。

图4  x = 7.5 m处波浪水平振荡速度幅值U_w垂向分布

Fig. 4  Vertical distributions of horizontal wave orbital velocity amplitude U_w at x = 7.5 m

与无植物工况相比，不同波浪条件下U_w在工况b和c的植物冠层内的垂向平均模拟值分别衰减了3 %和2 %（W1），6 %和8 %（W2），10 %和15 %（W3），9 %和17 %（W4）。植物的存在对波浪具有明显的阻碍作用，且植物引起的波浪速度衰减幅度与植物密度和波高均呈正相关关系。当波高较小时，水动力相对较弱，植物密度对速度衰减的影响较弱。反之，随着波高的增大，水动力增强，增大植物密度会引起更显著的速度衰减幅度。

垂向均匀密度植物影响下，如图4b、4c所示，植物的阻碍作用使得下层水体流速显著小于上层。在冠层顶部（z/h_v = 1），存在显著的流速梯度，出现U_w峰值，且该值随着植物密度和波高的增大而增大。为进一步分析该现象，图5绘制了工况a-W3和c-W3条件下x = 7.5 m处的波浪最大流速U_max和最小流速U_min的垂向分布情况。在工况c中，垂向上U_max在冠层顶部附近取得最大值，该值甚至大于无植物工况，说明植物的存在对该处的正向流速有一定的增幅。在冠顶处一个完整的波浪运动周期内，水分子的运动轨迹是在冠顶上方与波浪传播方向相同，为正向流速；而冠顶下方与波浪传播方向相反，为负向流速。冠顶处的植物拖曳力和水体惯性力会引起较大的雷诺剪切应力^［18］，增大了冠顶上方的正向流速，表现出冠层顶部的U_max极值；而冠顶下方的负向流速被植物拖曳力减小，使得冠顶附近的水分子运动出现周期不对称的现象，表现出正向的平均流速。该现象可等效为斯托克斯漂移，与粗糙底床上波浪边界层的流速类似，这也使得很多数学模型中将植物冠层概化为用拖曳力和惯性力模拟的粗糙多孔介质。

图5  工况a-W3和c-W3条件下x = 7.5 m处波浪最大流速U_max和最小流速U_min垂向分布

Fig. 5  Vertical distributions of maximum and minimum wave velocity (U_max and U_min) at x = 7.5 m in test a-W3 and c-W3

在工况d垂向不均匀密度植物的影响下，植物密度由水体上层向底部逐层增加，导致水体流速随之逐层减小，且在每个冠层顶部附近均有较大的流速梯度（图4d）。同时，在工况d中，各冠顶上方均有更高植物的存在，此处正向雷诺剪切应力对流速的增强作用会被上方植物的拖曳力作用抵消，因此在图4d中，各冠顶附近没有出现流速极值。

2.1.2　紊动能分析

图6为断面x = 7.5 m处紊动能k_T（$=\frac{\mathrm{1}}{\mathrm{2}}\left(\overline{{u^{\text{'}}}^{\mathrm{2}}}+\overline{{v^{\text{'}}}^{\mathrm{2}}}+\overline{{w^{\text{'}}}^{\mathrm{2}}}\right)$，其中$u^{\text{'}}$，$v^{\text{'}}$，$w^{\text{'}}$分别为x，y，z向的紊动流速）的垂向分布。无植物（工况a）模拟中，R_MSE较小（< 0.35 cm²·s^-2），但由于紊动强度较弱，模拟值与实测值的差距相对较大（I_A > 0.62）。与之相比，植物工况b~工况d模拟中，对波高较小的W1和W2，模拟精度较高（如工况c-W2中，R_MSE = 0.87 cm²·s^-2，I_A = 0.96）；而对波高较大的W3和W4，模拟精度较低（如工况d-W4中，R_MSE = 3.06 cm²·s^-2，I_A = 0.42），且模拟值小于实测值。这可能是由于波高较大、水动力作用较强时，物理模型试验中的植物会随波浪的传播发生轻微的晃动；而在数值模拟中的植物被设定为固定不动的理想圆柱体，未考虑其发生的晃动。

图6  x = 7.5 m处紊动能k_T垂向分布

Fig. 6  Vertical distributions of turbulent kinetic energy k_T at x = 7.5 m

图6a中，k_T由上层水体至下层逐渐减小。无植物影响时，波浪传播过程中，水面附近产生相对较大的紊动，且紊动由水面至底面逐渐减弱。与无植物条件相比，植物的存在使得水体的紊动明显增强，表现为工况b~工况d的k_T显著大于工况a。

如图6b、6c所示，在统一高度植物影响下，尾流紊动在植物的中部偏上位置处取得最大值，然后向冠层顶部和底部减小。在冠层下方，所有波浪条件下，工况c的k_T模拟值均比工况b大30 %~70 %，说明尾流紊动强度与植物密度正相关。而在植物垂向密度发生变化的工况d中，一方面，由水面至底部，植物密度逐渐增加，尾流紊动也随之增强；另一方面，尾流紊动的最大值出现在不同高度植物的中上部。两种效果的叠加使得图6d中k_T模拟值的垂向梯度较小。

选取工况c-W3条件下两根植物之间的三个测点，分析一个波周期内k_T和流速u的相位变化。如图7d所示，植物B位于植物A正前方5 cm，测点1位于植物A正前方7 mm，测点2位于两根植物中点，测点3位于植物B正后方7 mm（植物A和B半径均为4 mm，规定波浪传播方向为正）。对比三个测点的纵向流速u（θ），测点2的流速幅值显著大于测点1和3。当水流流向植物时，在植物前缘由于其阻流作用水流流速迅速减小并向两侧分散开来，而后在植物后方汇聚，形成尾流区，在尾流区至下一植物前缘区域流速逐渐增大至稳定。由图7d可知，当波浪流速为正向时，测点1位于植物A的尾流区，测点3位于植物B的前缘；当波浪流速为负向时，测点1位于植物A的前缘，测点3位于植物B的尾流区。因此，测点1和3始终处于流速减小区域，其流速幅值小于测点2。

图7  工况c-W3条件下，一个波周期内波浪流速u、紊动能k_T在不同测点处的相位图和三个测点位置示意图（y= 0.125 m，z = 0.3 m，尺寸不按比例）

Fig. 7  Wave velocity u and turbulent kinetic energy k_T versus wave phase at different measurement locations in test c-W3, and schematic of three measurement locations（y= 0.125 m, z = 0.3 m, not to scale）

如图7a、7c所示，在测点1和3处，波周期内紊动能k_T（θ）出现两个峰值，且最大峰值显著大于第二峰值。最大峰值与纵向流速u（θ）的峰值/谷值之间存在相位差，而第二峰值基本与u（θ）的峰值/谷值同时出现。在植物前缘（波浪流速为正时测点3，为负时测点1），植物的阻流作用使得该区域出现较大的流速梯度，进而产生较大的紊动，表现为测点1的u（θ）谷值和测点3的u（θ）峰值对应的k_T（θ）第二峰值。而尾流区（波浪流速为正时测点1，为负时测点3）的紊动需要经历在植物周围形成完整且稳定的涡旋后再脱落至水体中的过程，因此尾流区的k_T（θ）峰值滞后于流速极值，表现为测点1中k_T（θ）最大峰值滞后于u（θ）峰值约30°，测点3中k_T（θ）最大峰值滞后于u（θ）谷值约50°。

2.2　泥沙悬浮分析

2.2.1　Flow–3D原始泥沙模块模拟结果

采用Flow–3D原始泥沙模块（见1.1.2节），模拟W3波浪条件下泥沙悬浮过程，结果如图8所示。工况a中，模拟的垂向泥沙浓度分布与实测值吻合较好（R_MSE = 0.019 kg·m^－3，I_A = 0.99），说明该模型可以较好地模拟出裸床（无植物）条件下的泥沙起动及悬浮过程。然而，植物影响下的模拟结果却与实测值存在极大差别。工况b~工况d的模拟值远小于实测值，三组的评估指标R_MSE > 0.21 kg·m^-3，I_A < 0.61。含植物水流数值模拟中，植物的存在减小了流速（图4），由式（7）和（8）可知，床面切应力、Shields数和泥沙挟带上升速度均减小，泥沙起动量降低，表现出图8中含植物工况泥沙浓度模拟值小于无植物工况。而在物理模型试验中，植物的存在显著增强了水体紊动（图6），进而引起大量的泥沙悬浮，使得含植物工况泥沙浓度实测值大于无植物工况。由此可见，以切应力进行底部泥沙起动及悬浮判断的Flow–3D原始泥沙模块无法精确模拟植物影响下的泥沙运动。

图8  W3条件下x = 7.5 m处泥沙浓度C垂向分布

Fig. 8  Vertical distributions of sediment concentration C at x = 7.5 m in test W3

2.2.2　改进的Flow–3D泥沙模块模拟结果

为精确模拟含植物床面泥沙悬浮过程，利用基于水体紊动能改进的泥沙模块（见1.2节），对各工况下的泥沙运动过程进行计算。在计算过程中，基于对本研究水槽试验数据的分析，对式（17）中参数δ在床面附近取值为0.63较为合理^［16］，该值与Zhang等^［10］的推荐值（0.76）相近；对式（17）中参数C_w取值为0.015，该值在Tang等^［15］的推荐范围内（0.02 ± 0.01）。

图9为泥沙改进模块的计算结果，各工况下模拟结果均较好（R_MSE < 0.1 kg·m^－3，I_A > 0.89），说明该模型能够合理预测波浪条件下含植物床面泥沙起动及悬浮过程。虽然流速U_w和紊动能k_T的模拟结果存在部分偏差（图4、图6），但泥沙模拟值与实测值吻合较好，I_A值相对较高。这主要有两方面原因：一方面，物理模型试验中，用三维声学多普勒测速仪（ADV）测量流速参数时（采样频率64 Hz）对环境噪声的敏感性较高，且易受植物晃动影响；而光学后向散射浊度传感器（OBS）测量泥沙浓度时（采样频率8 Hz）不易受环境噪声影响，精度较高。另一方面，泥沙输运的模拟效果很大程度取决于近床面泥沙的起动及悬浮，而物模试验中植物晃动主要影响上层水体的紊动，对近床面的紊动影响很微弱。因此，数值模拟中，在式（17）中选择合理的参数δ、C_D和C_w使得改进后的泥沙模型充分考虑到床面剪切紊动和植物尾流紊动后，可得到较好的泥沙输运模拟效果。

图9  x = 7.5 m处泥沙浓度C垂向分布

Fig. 9  Vertical distributions of sediment concentration C at x = 7.5 m

图9中，模拟含植物床面（工况b~工况d）的泥沙浓度显著大于无植物床面（工况a），且波高越大越明显，与实测数据一致，证明刚性植物引起的尾流紊动对泥沙起动及悬浮的促进作用十分显著。在近底层（z/h_v < 0.5），工况c的泥沙浓度显著大于工况b，主要是因为高密度植物引起的尾流紊动更强，使得更多的泥沙起动。而在水体上部（z/h_v > 0.5），工况b和c的泥沙浓度垂向分布结构相近，主要是因为两种工况下紊动能k_T垂向分布结构相似（图6b、6c）。在工况d中，底层水体的泥沙浓度略大于工况c，而上层水体则相反。虽然工况c和d植物密度相同，但两者的k_T垂向分布结构差异较大（图6c、6d），导致了不同的泥沙浓度垂向分布规律。

3 结语

本文基于Flow–3D构建了含刚性沉水植物的三维水槽模型，改进了泥沙计算模块，模拟了不同波高的波浪通过植物区时水动力特征及泥沙输运过程，利用均方根误差（R_MSE）和一致性指数（I_A）证明模拟结果与实测数据较为吻合。主要结论如下：

（1）植物的存在使得波浪流速减小，其衰减幅度与植物密度正相关。在垂向密度均匀的植物群中，冠层顶部正向流速被增大；而在垂向密度不均匀的植物群中，由于不同高度植物的相互影响，该现象并不明显。

（2）植物的存在显著增强了水体紊动，其增幅与植物密度正相关。植物尾流紊动在植物中上部取得最大值。对单根植物而言，在其植物前缘和后方尾流区均属于流速减小区。当波浪通过植物区时，在每根植物附近一个周期内水体紊动能存在两个峰值，最大峰值与尾流区的涡旋有关，且与流速的峰值/谷值之间存在相位差；第二峰值与植物前缘的阻流作用有关，且与流速的峰值/谷值同步。

（3）在Flow–3D原始泥沙模块中，仅利用床面切应力计算泥沙挟带上升速度。将其应用于含植物床面时，没有考虑到植物尾流紊动对泥沙悬浮的促进作用，因此模拟出的泥沙悬浮量显著偏小。本研究从紊动能的角度改进了泥沙模块，在模型中即考虑植物尾流紊动，也考虑床面剪切紊动。改进后对泥沙悬浮的模拟评价指标由R_MSE > 0.21 kg·m^－3，I_A < 0.61提升为R_MSE < 0.1 kg·m^－3，I_A > 0.89，显著提高了含植物床面泥沙悬浮的模拟精度。

作者贡献声明

陈 明：构思框架，构建模型，处理数据，撰写、修改文稿。

娄 厦：构思框架，模型指导，论文思路，文本修改、确定文稿。

刘曙光：指导研究方案和论文撰写，全文审阅。

RADNAEVA Larisa Dorzhievna：数值模型技术支持。

NIKITINA Elena：数值模型框架优化。

CHALOV Sergey Romanovich：数值模型技术支持。