岸坡植被对复式河道水动力特性影响的数值模拟

郝由之，假冬冬,2*，张幸农，陈长英，吴 磊，杨 俊

(1.南京水利科学研究院 港口航道泥沙工程交通行业重点实验室，江苏 南京 210029；2.长江保护与绿色发展研究院，江苏 南京 210098)

天然河流大多是由主槽和漫滩组成的复式河道，枯水期主槽承担泄流功能，洪水期水流溢出主槽，由主槽和漫滩共同泄流[1]。而在自然河流近岸临水区，即主槽岸坡，经常会观察到乔木和灌木等植被。传统观点认为，岸坡植被不仅可以减缓水流速度，减弱河岸冲刷侵蚀，避免崩岸频发，从而稳定河势，还能起到净化水体、美化环境的作用，是一种既经济又有效的生态治河措施[2]。然而在洪水期，岸坡植被在一定程度上会改变滩槽内部原有水流结构，同时植被和水流间的相互作用会使河道中部分水流的能量转化为植被附近产生的紊动能，从而降低河道行洪能力。因此，了解复式河道岸坡植被对水流的扰动影响，进而分析岸坡植被的利害关系，对于河道防洪、护岸工程设计及水生态保护与修复均具有重要的研究意义[3]。

关于植被水流研究，在水槽试验方面，Nepf[4]利用示踪剂显示了圆形非淹没植被阵列内的紊流形态变化。López[5]、Nezu[6]和闫静[7]等以木棒或条状薄片模拟淹没刚性植被，研究了水流的紊动特性。赵芳[8]和杨克君[9]等研究了不同植被类型条件下河道的水流特性。随着数值计算理论和计算机技术的发展，解析方法和数值模拟技术得到成功的应用，且与试验相比，易得到整个计算域的水流特性。Huai等[10]提出了预测平均流速垂向分布的分析模型，很好地模拟了含双层刚性植被河道水流的垂向分布。Wu[11]和Zhang[12]等在动量方程和k-ε模型中考虑水流对植被的阻力影响，建立了植被水流2维模型。但是，含植被的河道水流具有复杂的3维各向异性，因此，3维模型的建立非常必要。吴梦瑶等[13]将多孔介质理论应用于自由表面水流3维模型，较好地反映了有限植被群的尾流区流动特征。刘诚[14]将Wu[11]等将模型扩展到非正交曲线坐标系下建立了3维模型。以上3维模拟方法可实现植被对水流的宏观效应模拟，但对植被群内紊流特性的捕捉仍存在一定困难。槐文信[15]和Nadaoka[16]等采用大涡模拟技术建立了植被水流模型，较好地模拟了植被尾流和植被群内细部水流结构和大涡的发展过程，但当计算域较大，而细部涡流复杂区域占比相对较小时，采用大涡模拟运算量较大。

复式河道受自身断面形态特点影响，水流特性复杂，再加上河道内植被与水体间相互作用，水体流动更加复杂。然而，当前研究大多以河床或漫滩植被为研究对象，如：米云彤[17]和任俊韬[18]等考虑了河床植被对水流特性的影响；蒋北寒[19]和袁素勤[20]等探究了漫滩植被对流速分布的影响。而实际河流中，植被在复式河道岸坡广泛分布，其对河道行洪的影响尚不明确。以往研究也多以单个植被或植被斑块为研究对象，而对顺水流方向植被沿程覆盖的研究较少，且现有研究结果主要集中于植被对河道流速分布和紊流特性的影响，对二次流结构、河道流量分配及河床剪应力分布影响的研究尚不多见，而同时考虑这些影响能够更加全面认识植被对水流的扰动。因此，岸坡植被对复式河道水动力特性的影响还需进一步研究。

本文将植被作为固壁边界，并采用雷诺应力模型封闭紊流控制方程，建立3维植被河道水动力数学模型，相比于2维模型，可以获取河道横断面的水流特性；同时，与传统研究在动量方程中添加植被阻力项相比，本文方法可有效捕捉植被区域内部的紊流结构，同时又避免了使用大涡模拟时的高运算量。为了验证模型的准确性，对复式河道的水流特性进行了模拟，并将结果与Tominaga和Nezu[21]的试验数据进行比较。随后，利用该模型模拟了复式河道岸坡种植不同淹没状态的植被和不同河岸坡比条件下的河道水流结构，以揭示岸坡植被对复式河道水流扰动的影响。

1 方 法

1.1 数学模型构建

1.1.1 控制方程

利用CFD计算流体力学软件建立基于有限体积法的3维植被河道水动力数学模型，其中紊流控制方程采用连续性方程和雷诺时均纳维-斯托克斯方程[8]，并通过雷诺应力模型对紊流模型进行封闭。

对于不可压缩均匀流体，连续性方程与雷诺时均纳维-斯托克斯方程形式如下：

式（1）～（2）中：i,j= 1, 2, 3为空间指数，x1、x2、x3分别代表x、y、z方向，即横向、纵向、垂向；ρ为水的密度，取998.2 kg/m3；μ为水的动力黏滞系数，取1.002×10-3Pa·s；p为时均压力；Ui（Uj）和ui（uj）分别为在xi（xj）方向上的时均流速和瞬时流速；为雷诺应力项；gi为重力加速项。

紊流封闭采用雷诺应力模型，该模型通过求解雷诺应力输运方程对方程组进行封闭[22]。雷诺应力输运方程具体形式如下：

其中：

对流项为：

式中：xk为垂直于xi和xj的方向；uk为xk方向的瞬时流速。

紊动扩散项可通过Daly等[23]给出的广义梯度扩散模型进行计算：

分子黏性扩散项为：

剪应力产生项为：

压力应变项为：

式中，p′为单位时间内的压力。

黏性耗散项为：

雷诺应力输运方程中，i、j取1、2、3时，得到3个正应力方程，将3个方程相加得到紊动能的准确方程；而方程式中除了剪应项外，均含有2阶或3阶相关矩，因此，需给定边界条件使雷诺应力输运方程封闭，才能对紊动能、紊动耗散率等参数进行求解。

1.1.2 边界条件

模型计算域及边界条件如图1所示，边界条件类型主要包含速度进口边界、自由流出口边界、对称边界和壁面边界。复式河道中，水面线几乎与底坡平行，自由水面高程的起伏变化很小，在水深较大的情况下可以用对称边界条件近似代替自由水面，因此，自由水面采用刚盖假定法做近似处理。模型进口处设定为水流速度进口边界条件，并给定平均流速Uave。出口处设为自由出流边界条件，即出口处沿水流流向的速度梯度和紊动强度梯度为零。由于河道形态沿主槽中心对称分布，为了减少计算工作量，计算区域选取整个断面的一半，将主槽中心纵剖面设为对称边界。河床底部、滩地边界及刚性植被均设为壁面边界，黏性流动中壁面默认为静止无滑移边界条件，即Ui= 0。由于Tominaga等[21]研究发现相对于断面形态与水力条件，壁面糙率对水流结构的影响十分微弱[24]，因此，该模型中未考虑壁面糙率的影响。

图1 模型计算域及边界条件示意图Fig.1 Schematic diagram of model calculation domain and boundary conditions

1.1.3 植被固壁边界处理方法及网格划分

网格划分方法和网格质量对模拟结果有很大的影响，因此合理的网格划分可以提高模型的计算精度和计算效率。针对本文模型的特点，要保证在靠近河床和植被边壁的流速大梯度变化区域网格有较细的间距，而在主槽和漫滩中心的低速梯度区域网格有更宽的间距。此外，主槽和漫滩区水平面上划分为均匀的笛卡尔网格块，而岸坡斜面上划分为向植被圆柱体拉伸的“O”型网格块。此划分方法能够在保证捕捉植被区细部水流结构的同时减小网格数量，节省计算资源。模拟时采用了不同分辨率的网格验证了模型的收敛性和模拟精度，进而选取适当的计算网格。以岸坡种植非淹没刚性植被为例，其3维几何模型不同剖面的网格划分如图2所示。

图2 几何模型在不同剖面的网格划分Fig.2 Grid division of the geometric model in different sections

图2中：近床面和岸坡植被区网格细密，主槽和漫滩中心网格较粗；同时，在x-z剖面和植被区水平剖面图中，x1∶x2、z1∶z2及y1∶y2应小于网格尺度比1∶5。

1.1.4 求解方法

采用基于有限体积法的计算流体力学软件建立含植被河道3维紊流数值模型，方程离散时，压力速度耦合方式选择SIMPLE算法，压力的插值格式采用体积力加权格式，动量等其余量的离散格式采用2阶迎风插值格式。当各计算参量残差小于0.000 01且监测某点的参量不随迭代次数变化时，可认为计算收敛，迭代终止。

1.2 模型验证

采用Tominaga和Nezu[21]矩形复式河道紊流试验中工况 S-2的结果对所建数学模型进行验证。该试验中，水槽长12.5 m，宽0.4 m，高0.4 m；水槽底板为铁板，边壁为有机玻璃。河道主槽水深hm为0.08 m，滩地水深hf为0.04 m，主槽宽度b为0.2 m，总宽度B为0.4 m，水力参数见表1。

表1 验证工况相关水力参数Tab.1 Related hydraulic parameters of verification conditions

选择河道中心横断面作为代表断面，绘制了二次流、主流速及河床剪切应力分布，如图3～5所示。图3～5中，横、纵坐标轴均使用主槽水深hm进行无量纲处理，主流速使用最大流速Umax进行无量纲处理，河床剪切应力 τ采用平均剪切应力 τ进行无量纲处理。

图3 二次流分布验证Fig.3 Secondary flow distribution verification

由图3可知，试验和模拟结果在复式河道滩槽交接处均存在一对方向相反的二次流漩涡，二次流漩涡的影响范围从滩槽交界处开始延伸至自由表面，水平方向的影响范围也一致且漩涡的位置、大小和强度十分接近，故拟合效果较好。

由图4可知，试验与模拟结果规律相近，只在自由水面处略有误差，这是由于数值计算时将自由水面作为对称面边界处理，无法较真实地模拟自由表面涡，导致近水面流速误差相对较大。但数值模拟与试验结果均显示，复式河道滩槽交界处的主流速等值线向主槽一侧的自由水面凸起，滩槽交界区主流速有明显减小趋势，说明数值模拟结果仍能较好地反映断面主流速分布的主要特点。

图4 主流速分布验证Fig.4 Main flow velocity distribution verification

由图5可知，复式河道河床剪切应力分布的数值模拟与试验结果基本一致，均在滩槽交界处剪应力迅速减小并达到最小值，并在交界处向主河道区和漫滩区两个方向过渡时剪切应力值迅速增加至最大值，然后再分别减小，直至边界处达到最小值。

图5 河床剪切应力分布验证Fig.5 River bed shear stress distribution verification

1.3 工况设置

模型设置为复式梯形河道，如图2所示，模拟类似洪水期漫滩水流河道，设置主槽宽度为1.5 m，漫滩宽度为2.0 m，岸坡高度为0.4 m，水深0.8 m，流向长度30.0 m，河道坡降为0.1%。岸坡顺水流方向沿程布设刚性植被，植被以平行方式排布，间距为0.14 m，植被直径为20 mm。模拟工况设置了不同的植被淹没方式及不同的河岸坡比（河岸横向宽度与高度之比）进行对比，所有模拟工况中的进口平均流速Uave均取0.5 m/s，具体工况设置见表2。

表2 数值模拟工况Tab.2 Numerical simulation conditions

2 结果与讨论

2.1 流速分布与流量分配

图6为岸坡种植刚性非淹没植被时（z=0.6 m）处的河道水平面流速分布。由图6可知，左侧进口处流场为非充分发展阶段，该区域流速分布与单个圆柱绕流相似[25]，也说明了模拟方法可以很好地捕捉植被区细部水流特点，具有一定的精度。由于本文考虑植被在岸坡沿程覆盖的情况，因此重点关注水流充分发展阶段的水流结构，故下文河道断面水流特性数据也均取自充分发展阶段的断面进行分析。

图6 水平面流速分布Fig.6 Velocity distribution on the horizontal surface

图7为8种工况下的复式河道横断面流速分布。

图7 断面主流速分布Fig.7 Distribution of main flow velocity in section

对比图7（a）～7（d）发现：岸坡植被使得岸坡区流速显著减小，且流速减小区域与植被分布相关；其中，刚性非淹没植被所占岸坡范围最大，因此对滩槽交界区流速的减小作用最强。

图7（e）～7（h）分别为1∶1.5和1∶2.0坡比条件下，岸坡有、无植被时的断面主流速分布，与图7（a）和7（b）对比可发现：无植被情况下，主流速在滩槽交界处凸起，这与Xiao等[26]对不同坡比下复式河道水流研究得出的结果相似，这是由于此处边界凸起，引起了二次流的产生；但当岸坡植被存在时，滩槽交互区流速仍显著减小。

岸坡植被改变了河道断面流速分布，进而影响流量在主槽、河岸和漫滩3个区域的分配，图8为工况1～4的断面流量分配占比。由图8可知，在来流量相同条件下，岸坡植被的存在大大削弱了岸坡区过流量，但主槽和漫滩区过流量增加，这在一定程度上增加了河道的过流负荷。其中，河岸种植非淹没刚性植被影响最显著，与岸坡无植被相比，岸坡区流量分配由9.21%缩减到2.07%，减小了80%，使得主槽区流量分配增加了6%，漫滩区流量分配增加了13%，这些结果可为确定河道行洪负荷提供依据。

图8 断面流量分配Fig.8 Cross-section flow distribution

2.2 二次流分布

图9给出8种工况下复式河道的二次流分布，横、纵坐标轴均使用主槽水深hm进行无量纲处理。

图9 断面二次流分布Fig.9 Secondary flow distribution in section

图9（a）～（d）显示：同一坡比条件下，岸坡无植被时，滩槽交界区存在一对方向相反的二次流涡旋。当岸坡种植非淹没刚性植被时，二次流涡旋数量增加，滩槽交界区有两对方向相反的二次流涡旋；当岸坡交替种植淹没与非淹没刚性植被时，滩槽交界区有3个相邻且方向相反的二次流涡旋；当岸坡种植淹没刚性植被时，滩槽交界区有一对方向相反的二次流涡旋，但相比无植被时，二次流涡旋水平影响范围增加。此外，根据二次流流速矢量图例可发现，岸坡有植被时的二次流约为无植被情况的3倍，说明岸坡植被可以增加二次流涡旋的数量、强度及影响范围，其中非淹没刚性植被影响最显著。

由图9（a）、（e）及（g）可见，当河岸坡比逐渐减小时，涡旋水平方向的影响范围进一步减小，二次流的形状也逐渐模糊，这是由于坡角变缓后，主槽与滩地之间存在一个流速过渡区，滩槽间横向动量交换减弱，二次流强度也随之减弱。与图9（b）、（f）及（h）对比同样可发现，岸坡植被能够促使二次流涡旋的产生，并扩大二次流的影响范围，增加滩槽间横向动量交换，但随着坡比的减小，植被的影响会随之减弱。

2.3 河床剪应力分布

图10给出了8种工况下，对平均剪应力 τ进行无量纲化后的河床剪应力分布。图10中，河床剪切应力τ可由CFD计算后直接得出，再将 τ沿床面平均即可求得平均剪应力 τ。

图10 河床剪应力分布Fig.10 River bed shear stress distribution

由图10可见：各工况下河床剪应力分布趋势基本一致，均在主槽和滩地中心接近均匀分布，主槽和漫滩向岸坡过渡时河床剪应力逐渐增大，进入河岸边坡区又发生突变减小；但不同工况下，主槽和滩地中心河床剪应力值及滩槽交界区河床剪应力的变化幅度有所不同。

对比工况1～4可以发现：同一坡比条件下，岸坡无植被时，主槽和漫滩河床剪应力最小，无量纲化剪应力值约为1；岸坡种植植被后，主槽和漫滩河床剪应力均有所增大，可增大至1.1～1.3。而对于岸坡区无植被时，河床剪应力最大，无量纲化剪应力值约为0.7～1.0，该结果与Xiao等[26]对坡比1∶1复式河道河床剪应力的模拟结果基本吻合；有植被时，河床剪应力大幅减小，甚至趋近于零。在主槽和漫滩向岸坡过渡区，岸坡种植非淹没刚性植被时，河床剪应力增加幅度最大，交替种植淹没与非淹没刚性植被时次之，种植淹没刚性植被时最小。这说明非淹没刚性植被能够引起强烈的滩槽动量交换，从而使得该处河床剪应力急剧增加。

对比工况1～2和5～8可以发现，在岸坡有、无植被情况下，随着坡比的减小，主槽、漫滩及滩槽交界区的河床剪应力均逐渐增大，说明随着坡比减小，二次流引起的动量交换强度变小，对边界剪切应力分布的影响减弱，河床剪应力分布向均匀化趋近。

2.4 雷诺应力分布

由图11（a）可知，岸坡无植被时，雷诺应力在岸坡区呈凸起分布。与图11（b）～（d）对比可以发现：岸坡种植植被后，滩槽交界区雷诺应力分布与无植被时有很大差异，雷诺应力分布不再呈凸起状态，而是与岸坡植被分布相关；岸坡主槽侧存在一个紧邻植被的雷诺应力值增大区，且为正值；岸坡漫滩侧存在一个紧邻植被的雷诺应力值减小区，且为负值。这说明岸坡植被的存在显著影响了雷诺应力分布的变化，使得临近岸坡的主槽和漫滩与岸坡之间的横向动量交换更加剧烈。

图11 断面雷诺应力分布Fig.11 Reynolds stress distribution of section

2.5 紊动能分布

图12为1∶1坡比条件下4种工况的复式河道横断面紊动能分布，其中，紊动能k由平均摩阻流速U*进行无量纲化处理。

图12 断面紊动能分布Fig.12 Distribution of section turbulent kinetic energy

由图12可知：无植被时，滩槽交界区紊动能呈凸起分布；而岸坡种植植被后，同雷诺应力分布相似，该区域紊动能不再呈凸起状态，而是与岸坡植被分布相关。岸坡种植非淹没刚性植被时，紊动能在岸坡植被区两侧显著增加，无量纲化的紊动能值可达到20；岸坡交替种植淹没与非淹没刚性植被时，紊动能也在岸坡植被区两侧显著增加，但幅度比前者略小；岸坡种植淹没刚性植被时，紊动能在岸坡植被区两侧及植被顶部显著增加，增加幅度也小于前者。这说明岸坡植被的存在能够增加紧邻植被周围区域的紊动交换，而植被区域内部紊动能会显著减小。米云彤等[17]对植被群周围水流特性进行测定，其植被群处断面紊动能结果与本文结果趋势相同。

3 结 论

本文采用雷诺应力模型封闭紊流模型，建立了3维植被河道水动力数学模型，在模型验证相似的基础上，利用所建模型研究了岸坡刚性植被对水流的影响，主要研究结论为：

1）岸坡植被增加滩槽间动量交换，使得岸坡区的流速、流量及河床剪应力减小，而二次流涡旋数量、强度和范围均有所增加；此外，随着河岸坡比减小，无植被河道滩槽间动量交换逐渐衰减，但在坡面种植非淹没刚性植被后，滩槽间动量交换的衰减程度有所减小。植被区域内紊动有所减弱，而临近植被的周边区域紊动交换剧烈，且岸坡种植非淹没刚性植被时影响最显著。

2）根据断面流量分配结果，可确定岸坡植被对行洪负荷的影响，而河床剪应力结果又可与河岸土体起动建立联系，进一步为坡脚冲刷及河岸稳定提供判据，因此该模拟方法和结果有助于评估和权衡植被的利弊影响，进而制定适宜的河岸带生态护岸植被配置方案，使得植被发挥到最大的生态效益和工程安全效益。