刚性植被密度对岸滩作用的数值模拟*

肖苡辀，张国志，冯先导，陈文峰

(1.中交第二航务工程局有限公司，湖北 武汉 430040；2. 长大桥梁建设施工技术交通行业重点实验室，湖北 武汉 430040；3. 交通运输行业交通基础设施智能制造技术研发中心，湖北 武汉 430040；4. 中交公路长大桥建设国家工程研究中心有限公司，湖北 武汉 430040)

岸滩修复工程通常会利用植被来增加岸坡的稳定性，为水生动物提供遮荫和避难所。水生植被常存在于河岸、海岸边，可显著改变局部流场、泥沙输移、河岸稳定性和水生栖息地[1-5]。了解植被对水流、泥沙输移的影响是保护港口海岸岸线的关键。

近年来，学者对植被层内外的水流悬移质泥沙输移规律及床面推移质泥沙起动特性进行了系统研究。张英豪等[6]通过试验研究了植物对水流和泥沙运动的影响、植物的消波与捕沙作用，认为合理地布置植被可以减少泥沙输运；Vastila等[7]通过试验发现清水中的泥沙沉降与芦苇植被水流中的泥沙沉降基本一致；薛万云等[8]对刚性植被环境中床面形态变化进行了试验研究，结果表明淹没植被条件下动床床面泥沙运动可分为静止、部分起动及全部起动3个阶段。

综上，植被对岸滩的稳定性分析主要采取试验的研究方法，缺乏数值模拟植被密度对促淤保滩作用方面的研究。本文选取圆柱代替刚性挺水植被，采用数值模拟方法探索植被水流中床面变化形态规律，分析植株周围水流流速、床面切应力与植被密度的关系，对促进植被-水流-泥沙研究的发展具有重要价值。

1 数值模拟方法

1.1 控制方程

假定水流为不可压缩黏性流体，则水流计算可运用Navier-Stokes方程描述[9]：

(1)

(2)

式中：ρ为流体密度；ui为局部时均速度分量；u为速度矢量；t为时间；μ为流体动力黏滞系数；p为压强；Fi为微元体上的体力。

本研究采用的泥沙冲刷模型基于文献[10]。Soulsby-Whitehouse方程[11]用于计算临界Shields数θcr，i：

(3)

(4)

式中：ρi为第i相泥沙密度；ρf为流体密度；di为颗粒直径；μf为流体动力黏滞系数；g为重力加速度。

局部的Shields参数θi是基于局部床面切应力τ计算的：

(5)

其中τ考虑床面粗糙度，采用壁面定律进行计算。

利用夹带系数(0.005)来衡量冲刷速率，并与试验数据进行拟合。沉降速度控制着Soulsby沉降方程。公式(6)～(8)用来计算河床单宽的体积输沙率：

(6)

Φi=βMPM，i(θi-θcr，i)1.5cb，i

(7)

(8)

式中：βVR，i、βMPM，i和βNie，i一般为0.053、8.0和12.0；cb，i为第i相的体积分数；Φi为无量纲推移质输沙率，假定为0.05。容积床载输送率qb，i表示为：

(9)

根据粗糙度修正后的壁面平均速度定律为：

(10)

式中：up是在某点j处的主流流速；k为von Karman常数，为0.418；u为与底部剪切应力有关的剪切速度；cμ为常数；kj为点j处的紊动能；ΔB为粗糙度函数；zj为点j的垂向高度。

1.2 流场设置与边界条件

本研究基于袁代亮[12]的试验设置建立模型，通过CFD软件计算流场，旨在验证模型的准确性。该验证模型长3.0 m、宽0.5 m，水深0.24 m，流速0.227 8 m/s，主要监测其垂向流速分布。在模型得到验证后进行不同植被密度下的水流-泥沙计算(图1)，具体工况参数见表1。计算区域原点设于水槽进口底板中心点，整体计算区域长3.0 m、宽0.5 m、高0.5 m。刚性植被采用圆柱代替，直径为0.1 m，第1排植被(植被圆心)位于距离水槽进口1 m，植被纵向间距为0.1 m，植被横向间距依据植被密度不同而不同(从工况1～4，横向间距分别为0.125、0.100、0.083、0.071 m)。泥沙中值粒径为0.006 m，密度为2 650 kg/m3，厚度为0.23 m，铺设在水槽床面上。

图1 计算模型

表1 试验工况设计

模型采用流量进口边界作为进口边界条件，压力出口边界作为出口边界条件。壁面边界用于水槽底板和边壁。渠道顶部空气入口设定为对称边界，即默认无流体穿过该边界。

1.3 计算方法

基于有限差分法采用RNGk-ε紊流模型来求解控制方程，对流项采用二阶迎风格式，扩散项采用二阶中心差分格式。采用TruVOF方法缩减收敛时间，该方法仅计算流体单元，不计算空气单元，且对自由液面的描述更为准确。通过FAVOR法(fractional area volume obstacle representation)对计算区域划分网格。当控制方程的残差小于1×10-5时，认为计算结果收敛。

2 结果分析

2.1 试验验证

为了更好地研究植被存在时水动力结构和悬沙运动特性，需要对数值模拟的可靠性进行验证。基于袁代亮的试验结果进行数值模拟，从而对比分析数值模拟的可靠性。本研究基于试验设置和结果进行数值模拟，结果见图2。图2为垂向流速分布模拟值与试验值对比，结果表明模拟值与试验值吻合度较好，该数值模拟方法能够有效模拟植被水流工况。

注：z为垂向高度，H为该工况水深。图2 垂向流速分布模拟值与试验值对比

2.2 流速分布

基于速度场的平均水力特性研究植被对水流结构的影响。图3为各工况水平截面z/D=15的速度分布。由图3可知，植被的存在阻碍了水流的流动，而植被间的流速则增大。值得注意的是，从工况1～3，随着植被密度的增加，流速逐渐增大；但比较工况3和4，植被密度的增加反而降低了流速。这个现象是在先前研究中出现过的[13]，其原因是工况1～3中，卡门涡街在动量运输中占有主导作用；而在工况3和4中，卡门涡街已经被密度较高的植被耗散了，以至于涡街并未达到下一排植被。

图3 各工况水平截面zD=15的速度分布

2.3 床面切应力

床面切应力影响着泥沙输移，因此在工程中具有重要意义。床面切应力τ计算公式如下：

(11)

式中：Ud是水深平均流速；Cf是表面摩擦系数。

图4呈现了各工况横断面和工况4沿程横断面床面切应力分布情况，其中图4a)为工况4的监测断面分布，由虚线表示。由图4可以看到，植被水流中的床面切应力分布呈现波浪形状，这说明植被的存在降低了平均流速，消耗了水流动量，增加了阻力，从而降低床面切应力。由图4b)可知，从工况1～4，床面切应力随着植被密度的增加而出现先增大再减小的趋势，其中工况3的床面切应力最大。由此可知，一味地增大植被密度并不一定能够保护好滩地，而是要到一定密度之后才会降低床面切应力从而保护床面。

由图4c)可知，工况4断面的床面切应力呈现沿程逐渐增加的规律。在植被群上游断面(断面x/D=93)，床面切应力分布较为均匀；当水流经过第一排植被后，床面切应力开始呈现波浪形状分布；随后在植被群内，床面切应力沿程先增大后减小；直至水流经过所有植被后，床面切应力逐渐减小。比较图4c)的横断面x/D=93和x/D=175床面切应力沿程变化规律显示，植被的存在减小了床面大部分的切应力，具有促淤保滩的效果。值得一提的是，在本研究中，横断面靠边壁处的床面切应力经过植被群之后是增加的，这是由于计算区域的边界受到限制而导致的。

图4 床面切应力分布对比

2.4 床面形态变化

前文已对植被对水流结构影响进行了定量分析，现对床面泥沙的影响进行直观展现。图5展示了各工况床面形态变化，较为直接地呈现了植被水流对床面泥沙的作用与影响。由图5可见，在单根植被周围形成了冲刷孔，在单根植被下游形成沙脊。这是由于水流在遇到植被时形成圆柱绕流，其生成的马蹄涡向下淘刷形成冲刷孔；单根植被下游近区存在滞水区，滞水区流速较低，泥沙容易在此沉积，所以形成了沙脊。同时，在图5中发现，从工况1至工况3，随着植被密度的增加，植被区域形成沙脊逐渐明显；而从工况3到工况4时，植被密度的增加使整体形成沙脊的现象不那么明显，而是以单根植被下游近区的局域沙堆出现。这说明了涡街被较高密度的植被耗散了，以至于涡街并未达到下一排植被，只能在单根植被下游近区局部形成沙堆，从而减少了植被群下游床面泥沙的冲刷。

图5 各工况冲刷情况

3 结论

1)与试验结果对比表明，过流流速与试验结果较为吻合，该数值模拟方法可以有效地进行水流流场复现。

2)随着植被密度的增加，流速呈现先增大后减小的规律，其中原因是由于卡门涡街在较高密度的植被中容易耗散，以至于涡街并未达到下一排植被。

3)床面切应力沿程变化规律表明，植被的存在减小了床面大部分的切应力，具有促淤保滩的效果；随着植被密度的增加，床面切应力呈现先增大后减小的规律。

4)床面形态变化规律显示，在单根植被周围形成了冲刷孔，在植被下游形成沙脊；当植被密度增加到一定程度时，整体形成沙脊的现象不那么明显，而是以单根植被下游近区的局域沙堆出现。