基于有限体积法的不同植被茎粗下坡面流水力特性数值模拟

徐雪峰，张升堂，张景洲，周建森

（山东科技大学地球科学与工程学院，山东青岛266590）

0 引言

近年来，有学者对前人的研究进行分析总结，最终将坡面流定义为以降雨为主导，重力为辅助动力，流经不同坡面产生不同运动状态的浅层水流［1，2］。对于现实自然环境而言，水力侵蚀是土壤坡面侵蚀中最为广泛的自然现象（除风力侵蚀以外）。坡面水力侵蚀的演化与水流强度、坡面土壤的抗剪切能力两者的非协调发展机制有关，而坡面植被的存在，使得坡面水流形态更加特殊与复杂化，一方面使坡面水流内部结构发生改变，同时植被根系具有固结土壤作用，增强了土壤抗蚀性［3，4］，从而成为根治坡面水土流失的重要生物措施。

在国内外学者的研究工作中，坡面流植被区域的流动特性和湍流特性受到了广泛的关注。大量文献提及植被分布模式对坡面流阻力的影响。Wang等［5］重点研究了植被、水流和泥沙输移间的相互作用。植被密度和布局对流速有不同程度的影响，流阻随着植物密度的增加而增大。Zhao 等［6］认为坡面植株的根茎能有效的降低泥沙的输送能力，减少水土流失。张冠华等［7］通过试验研究植被的分布格局对坡面水流阻力的影响，对不同格局的阻力系数进行对比寻找规律。Zhang 等［8］为研究植被对水流的影响，通过理论分析和实验验证，得出Manning糙率系数n随雷诺数Re的增加而增大，且当植被覆盖度越大时，Manning 糙率系数n随雷诺数的变化率越大。姬昌辉等［9］采用水槽试验研究了在不同水流及种植间距下的淹没植被曼宁糙率系数的变化特征，结果表明糙率系数与植被间距呈负相关性。Ding 等［10］通过试验得出植被分布在试验段的下部比分布在试验段中部及上部的拦沙效果更好，证明了植被的分布模式对明渠水流的水力特征有着重要的影响。但大都侧重于实验模拟，不能准确得到植被周围的水流流动结构。

近年来也有一些研究人员对不同排布方式的植被进行了数值模拟的研究。Zhao和Huai［11］利用大涡模拟（LES）模型研究了明渠中不连续和淹没植被斑块对水流湍流的影响。Lima等［12］和Weber等［13］采用计算流体动力学（CFD）技术研究了明渠水生植被周围的水流流动结构。Ghani 等［14］采用数值模拟的方法，研究了圆形和交错植被斑块周围的流动情况，利用不同植被密度和流速的实验数据对数值模型进行了验证。Anjum等［15］运用CFD 工具FLUENT，对一个纵向不连续和垂直两层植被占据河道半宽度的水流进行了数值模拟。结果表明林隙内的流速明显慢于植被斑块内的流速，局部和不连续植被对水流结构和阻力有较大的影响。占据半宽度的垂直分层和不连续植被的存在显著地影响了纵向和横向的水流。Anjum和Tanaka［16］采用k-epsilon 模型，通过覆盖整个河道宽度的间断垂直分层植被来阐明河道水流的流动结构。然而这些研究大多是在相同粗度植被条件下的水流，但实际坡面植被复杂多样，茎秆粗细程度不同，即使同一坡面上的植被，在不同的生长阶段其茎秆粗度是不断变化的，而且现实中洪水水流和灌溉水流在植株间的流动受植株直径影响极大。

因此，对坡面植被植株粗度进行系统性的研究具有重要意义，有助于通过模型尺度对实际含植被坡面的水流阻力进行相关的预测，从而达到提高实际坡面的阻水能力，降低坡面水流的冲刷能力以及减小水力侵蚀的目的，为坡面水流阻力的研究提供价值性参考。为分析坡面植被茎杆粗度的变化对坡面流水力特性的影响，通过试验分析研究不同流速情况下的水流阻力及流态变化规律，从而得出植被的茎杆粗度的变化对坡面流水力特性影响的定性及定量研究论点。本文的研究目的是建立一个数值模型，数值模拟是由CFD 工具FLUENT 通过采用Standardk-ε紊流模型进行的，通过均匀分布的非淹没植被，分析坡面流的平均速度分布和湍流特性。研究还重点分析了植被间的水流分布特点。

对圆柱状模拟植被在坡面径流冲刷试验和数值模拟相结合的基础上，旨在研究和分析：

（1）数值模拟方法对于含植被坡面径流冲刷试验是否具有可行性；

（2）不同粗度变化的植被对流速的影响；

（3）研究植被间特定位置和特定截面的局部流动特性。

本数值模拟试验以统一的植被行列走向与水流方向所成夹角θ= 90°、相邻植株的行列间距a×a=60 mm×60 mm、高度为0.15 m的圆柱状塑料棒均匀插入相应钻孔口径的试验底板。为更加直观的反映植被茎杆粗度对坡面水动力特性的影响，本数值模拟试验在坡度s为0 的条件下，设置3 组茎杆粗度d（3、4、5 mm）的植被，16 个单位流量q进行了模拟计算，基于有限体积法，采用Standardk-ε紊流模型，建立数值分析模型，模拟刚性植被粗度变化的条件下坡面流水流阻力的影响，对三种粗度的刚性淹没植被明渠水流特性进行了研究。

1 数学模型

1.1 控制方程

在重力驱动下，有自由液面的坡面流动中，连续性方程为：

动量守恒方程：

式中：ui表示流速u在xi方向上的分量；ρ为流体密度；p为平均压强；ν为运动黏度。

1.2 湍流模型

采用Standardk-ε模型［17］，不考虑流体压缩时，Standardk-ε湍流模型的湍动能k和湍流耗散率ε的约束方程如下：

k方程：

ε方程：

式中：xj为坐标分量；ui、uj为平均相对速度分量；Pk是湍动能生成项，定义为为湍流耗散率，且ε=为动力黏滞系数；涡旋黏性系数Cμ= 0.0845。σk和σε分别为k和ε的普朗特数，σk= 1.0，σε=1.3；模化常数C1ε= 1.44，C2ε= 1.92，C3ε= 1.0。

2 三维湍流数值模拟分析

2.1 物理模型

本数值模拟试验模型参照实验室试验仪器建模，实体试验仪器为长5 m，宽0.4 m，高0.3 m 的循环矩形水槽，水槽共分三个区域：上游稳水区、试验铺设区、下游量水区，其中试验铺设区长度为3 m，上游稳水区及下游量水区各为1 m，水槽共分为五个断面，每个断面之间的间距均为0.75 m。在试验铺设区铺设有机树脂板，在板中钻孔插嵌圆柱形塑料棒模拟植被区。数值模拟试验可以通过初始化设置达到计算所需要的水流参数条件，因此不需要单独划分上游稳水区和下游量水区，所建物理模型只包含实验铺设区。

2.2 网格划分及网格无关性验证

利用ANSYS FLUENT-Mesh 绘制计算网格，网格划分情况如图1所示。网格数量达到200 多万个，正交质量平均值达到0.9 以上，单元质量平均值达到了0.7 以上，表明网格质量较好，满足数值模拟运算要求。对数值计算的网格进行了无关性验证（以4 mm 径粗，0.08 m 水深为例），结果如图2所示。随着网格数量的增加，水流平均流速不断增加，当网格数达到2 365 896 时，进一步增加网格数量，水流平均流速的增量小于0.5%，可认为此时的网格数已不显著影响数值计算的结果。综合考虑计算精度与计算效率，选用全流道网格总数为2 365 896的网格进行计算。对数值模型划分不同数量的网格数并通过计算所得水流平均流速的变化来检查网格数的无关性。当网格数对水流平均流速的变化影响小于1%时，认为网格数达到无关性要求。

图1 网格划分图Fig.1 Meshing diagram

图2 不同网格数下的水流平均流速Fig.2 Average flow velocity at different grid numbers

2.3 边界条件设置

①入口边界条件设置为质量入口，湍流强度选择5%，为一般水流的湍流强度；②出口边界为压力出口，出口压力设置为0，温度为300 k；③上表面利用滑移边壁来代替水气交界面，植被表面及壁面均采用标准壁面函数（这是更好地预测边界层中流动的必要条件），为无滑移边界，壁面边界条件设置为固壁边界。仿真和后处理在计算流体动力学（CFD）工具FLUENT 中进行。计算区域如图1所示。采用SIMPLE 算法来实现压力—速度耦合。在运算过程中，当所有残差都低于1×10-5时，认为计算结果是收敛的。在进一步的迭代中，残差没有变化，并且还检查了具有主流速的解，在后续的迭代中，该解不再明显变化。因此，使用这些标准，可以假定解决方案已达到稳定状态。

3 数值模拟结果分析

3.1 3 种植被茎粗在不同流量条件下粗度对流速的影响

该数值模拟依据图1所示模型及尺寸进行物理建模，试验数据来自Zhang［18］等人的试验资料，根据定床阻力试验要求，对坡面流水力学特性进行计算机数值模拟。该模拟研究了与之前试验数据相同流量条件下，不同植被茎粗情况下坡面流的情况下的流动状况。模拟结果通过FLUENT后处理软件对其进行整理分析，并与试验数据进行对比。不同植被茎杆粗度下试验数据与数值模拟数据水深h和平均流速v（h～v）关系曲线图如图3示。

图3 3种不同植被茎秆粗度h～v关系曲线图Fig.3 Three different vegetation stalk thickness h～v curve diagram

根据数值模拟结果发现，3种不同植被茎秆粗度的h～v关系总体呈现出随着流量的不断增大，水深随之增加，流速不断增大；当水深h较小时，流速随水深的变化较大，且不同植被茎秆粗度的速度差异较大，随着水深的增加，流速呈缓慢变化趋势。在同一水深h下，不同植被茎秆粗度情况下，水流流速不同；植被茎秆粗度越小对应的流速值也越大，3 种植被茎秆粗度的水流流速关系为：v3mm＞v4mm＞v5mm，这一规律说明坡面植被茎秆粗度的微小变化对水流流速有一定的影响，植被茎秆粗度越大，水流流速越小。

试验结果与数值模拟结果进行对比如图3 示，数值模拟结果比试验实测数据偏大，但保持在10%以内，因此可以认为数值模拟在坡面流水流阻力试验研究中的应用是可行的。对于坡面流而言，影响阻力系数的因素众多，在本次对比验证试验中，造成这种差异的原因可能是：在数值模拟试验中，雨滴打击力在试验中未涉及，因而可忽略不计；数值模拟试验中重力和水压力也已通过参数设置完成，水槽壁面条件设置为光滑壁面，为理想环境。从以上分析可见，误差产生的主要原因可确定为坡面流表面张力的作用以及光滑的壁面条件。

3.2 3 种植被茎粗在相同流量条件下水流时均流速分布

本数值模拟选取3、4、5 mm 3 种植被茎粗分为3 组，计算了在相同流量和水深植被淹没状态下的水流流速。并且每组模拟中各选取了16个测点。

测点分布如图4（b）所示，测点a～e取在模型中心线上，相距15 mm，其中a 和d 在植被后，b 和e在两植被中间；测点f～k 取在中心线沿y轴负方向平移15 mm；测点l～p 取在中心线沿y轴负方向平移30 mm。水流速度坐标为x轴，与速度垂直方向为y轴，纵坐标为z轴。每个测点取50 组数据，总共采集了2 400 个数据，保证了垂向速度分布的准确性。

图4 计算域示意图（单位：mm）Fig.4 Schematic diagram of calculation domain

3.2.1 纵向流速分布

图5～7所示分别为3、4、5 mm 直径的植被在相同流量条件下各测点水流流速的垂向分布。总体看来，处在两行植被中间的测点l～p纵向流速相对较大，靠近植被的测点f～k纵向流速次之，与植被同行分布的测点a～e纵向流速最小。

图5 3 mm直径植被测点纵向流速的垂向分布Fig.5 Vertical distribution of longitudinal velocity at 3 mm diameter vegetation survey point

图6 4 mm直径植被测点纵向流速的垂向分布Fig.6 Vertical distribution of longitudinal velocity at 4 mm diameter vegetation survey point

以（图5）3 mm 直径植被测点纵向流速的垂向分布为例。可以看出，两植被之间的测点a～e中，位于两植被中间位置的测点b 流速最大，测点d 位于植被正下游且紧邻植被，受到植被阻力影响最大，其纵向流速变化最小，相对于其他测点的纵向流速小的多。这是因为水流经过植被时做圆柱绕流，水流流经植被前受到植被的阻挡，迎流面压强势能逐渐增大，动能逐渐变小，因而速度逐渐变小。测点a 纵向流速次之，测点c 和e 速度相近，这是因为测点c 和e 均位于两植被之间，离植被较远，水流受到植被阻力影响较小，且由于两测点相距较近，可以忽略流速损失。处在两行植被中间的测点l～p中，测点n纵向流速最大，达到0.298 4 m/s，且为所有测点流速中最大。通过（图8）测点i、nyx剖面流速分布云图也可以看到这一现象，测点n 位于两株植被中间，流速表现最大。测点i所处位置更接近植被，流速表现小于测点n。我们还可以看出，从两株植被中间位置（测点n）到植被位置，流速表现为逐渐变小。这与植被的挡水作用有关。

图7 5 mm直径植被测点纵向流速的垂向分布Fig.7 Vertical distribution of longitudinal velocity at 5 mm diameter vegetation survey point

图8 测点（a～p）yx面流速分布云图Fig.8 Measurement point（a～p）yx surface velocity distribution cloud map

水流纵向流速随着到植被距离的增加先增加后减小。且随着植被粗度的增加，水流纵向流速逐渐变小。

3.2.2 垂向流速分布

考虑到植被对水流流动结构的影响较大，有必要对含植被水流的垂向流速变化进行研究。从3种粗度植被测点纵向流速的垂向分布（图5～7）可以很直观的看出，在模拟的植被区域内，速度的垂直变化是单调的，最小速度发生在河床。这与Zeng C［19］的测量结果相同。与Neumeier［20］的测量结果不同，Neumeier 的测量结果表明，最小流速发生在地层上方约5 cm处。主要原因是Neumeir实验中采用了变截面自然植被。最大植被密度出现在床上5 cm 左右，对应最小流速水平。位于不同xz剖面的测点a～e、f～k、l～p垂向流速分布呈现很明显的规律，垂向流速v（测点l～p）＞垂向流速v（测点f～k）＞垂向流速v（测点a～e），且流速分布差异非常明显。以3 mm 粗度植被纵向流速的垂向分布为例，垂向流速v（测点f～k）与垂向流速v（测点a～e）的差值达到0.150 7 m/s，增大168%，可以看出，在y方向距植被越远，垂向流速的增长率逐渐减小。通过（图10）测点a～pxz剖面垂向流速分布云图可以直观地体现出流速变化。

以（图5）3 mm 直径植被测点纵向流速的垂向分布为例，测点f～k 垂向流速分布几乎一致，植被淹没高度h=0.02 m 左右时出现拐点，垂向流速趋近于最大值（0.219 8 ～0.239 9 m/s）。测点l～p 的垂向流速分布一致性更强，植被淹没高度h=0.03 m 左右时出现拐点，纵向流速趋近于最大值（0.287 3 ～0.298 4 m/s），它们的流速垂向分布曲线近乎重合。这一现象通过（图10）测点a-pxz剖面垂向流速分布云图能够直观体现。由图10（b）和（c）可以发现f～k、l～p 测点所在xz剖面的垂向速度分布较为均匀，所以各个测点的流速垂向分布相差较小。与之相比，测点a～e 所在xz剖面的垂向速度分布变化较大，植被淹没高度h=0.02 m 左右时出现拐点，垂向流速趋近于最大值（0.089 5 ～0.141 8 m/s），如图10（a）所示，因此各个测点的流速垂向分布相差较大。

图10 测点a～p xz剖面垂向流速分布云图Fig.10 Vertical velocity distribution cloud diagram of the a～p xz profile at the measuring point

整体看来，与植被同行分布的测点a～e 垂向流速变化幅度最小（其中测点d 变化幅度最大，为0～0.147 9 m/s）。靠近植被的测点f～k 垂向流速变化幅度大于测点a～e（其中测点i 变化幅度最大，为0 ～0.239 9 m/s）。处在两行植被中间的测点l～p 的垂向流速变化幅度最大（其中测点n 变化幅度最大，为0～0.298 4 m/s）。通过（图9）测点a～pyz剖面垂向流速分布云图可以看出这种规律。且3组粗度植被的测点纵向流速的垂向分布图都呈现这种规律。通过进一步分析发现，这与植被对水流的紊动作用有关［21］。测点a～e 与植被同行分布，由于植被与水流交界面处水流紊动较为剧烈，紊动强度明显增加，其紊流强度垂向分布，在植物顶部的水流剧烈运动产生较大的能量损失是形成水流阻力的重要原因，因此测点a～e 垂向流速变化幅度最小。而测点f～k 与植被排布处较远，植被对水流的紊动作用减小，因此测点f～k 垂向流速变化幅度大于测点a～e。处在两行植被中间的测点l～p 离植被分布区最远，植被对水流的紊动作用最小，因而垂向流速变化幅度最大。

图9 测点（a～p）yz剖面垂向流速分布云图（单位：m/s）Fig.9 Measurement point（a～p）yz profile vertical velocity distribution cloud map

另外需要注意的是，在数值模拟数据结果中，在水槽底部有速度突增，Liu等［22］，Huai等［23］和Pang等［24］在植被水流的试验中也发现了同样的现象，Liu 等［22］认为水槽底部的速度突增可能是由植被底部的马蹄形涡旋或是结合涡旋引起的，马蹄形涡旋会使得水流从周围区域较快的绕过植被底部，从而导致植被底部附近流速突增。

4 结论

采用标准紊流模型，建立数值分析模型，进行三维数值计算，通过在流场入口设置为质量入口的方式，对矩形明渠中不同茎径连续植被水流的流动进行了模拟。计算所得的流速值与实测值吻合良好，验证了数值模拟方法的正确性和可行性。在建立数值分析模型的基础上，对坡面水流的流动过程中的水动力学特性的变化规律进行分析，可以得出以下结论。

（1）试验结果与数值模拟结果进行对比分析发现，模拟结果比试验实测数据偏大，但误差保持在10%以内，与试验值吻合良好，验证了数值模拟在坡面流水流阻力试验研究中的可行性。

（2）根据数值模拟结果发现，3 种不同植被茎秆粗度的h～v关系总体呈现出随着流量的不断增大，水深随之增加，流速不断增大；植被茎杆粗度增大，流速减小。这是由于植被茎杆粗度增大，相对应的总体覆盖度增大，水流与植被接触面增大，从而使水流阻力增大

（3）纵向流速分布规律：总体看来，处在两行植被中间的测点l～p 纵向流速相对较大，靠近植被的测点f-k 纵向流速次之，与植被同行分布的测点a～e纵向流速最小。水流纵向流速随着到植被距离的增加先增加后减小。且随着植被粗度的增加，水流纵向流速逐渐变小。这表明植被根部附近有利于沉积物的积聚，对植物群落以及水生生物提供营养物质环境和空气。

（4）垂向流速分布规律：在y方向距植被越远，垂向流速的增长率逐渐减小。从3 种粗度植被垂向分布图5～7 可以很直观的看出，位于不同xz剖面的测点a～e、f～k、l～p 垂向流速分布呈现很明显的规律，垂向流速v（测点l～p）＞垂向流速v（测点f～k）＞垂向流速v（测点a～e），且流速分布差异非常明显。测点a～e垂向流速变化幅度最小，测点f-k 垂向流速变化幅度大于测点a～e。处在两行植被中间的测点l～p 离植被分布区最远，植被对水流的紊动作用最小，因而垂向流速变化幅度最大。植被间隙的低流速有利于泥沙沉积物沉积，有助于坡面稳定，维持生态环境。□