水-草-土耦合边坡流场特性研究

陈 彬，姚 锡 伟，毛 朝 勋

(内江职业技术学院，四川 内江 641000)

0 引 言

随着绿色生态理念的深入人心，明渠或边坡中的植被流引起了大量研究者的关注[1-3]。当地表水流经植被区域时，植被按水流条件主要分为淹没和非淹没两种类型[4]。宋滢汀等[5]通过玻璃棒模拟刚性植被，采用三维激光测速仪对含植被明渠的水流流速和紊动强度进行测量，得到了含刚性植被明渠水流特性。陈甜甜等[6]将植被水流分为树干层、树冠层和自由水层，在Boussinesq涡流黏度理论和经典对数律流速公式的基础上，推求出淹没植被垂向流速分布解析解，并设计相应试验验证解析解的合理性。韩宇等[7]通过随机森林RF、BP神经网络、支持向量回归SVR和RBF神经网络等机器学习算法研究了明渠出露刚性植被流速分布，发现模拟结果与试验结果吻合良好。任姗等[8]利用人造草模拟柔性植被，通过PIV技术研究了柔性淹没植被的垂向流速、应力和紊流强度随曼宁系数的变化特性，发现应力的峰值出现在植被顶部，并与植被密度呈正相关。此外，一些研究者也通过数值模拟的方式研究含植被明渠或边坡中水流特性[9-11]。Gac[12]采用大尺度格子玻尔兹曼方法分别模拟了刚性和柔性植被下明渠垂向流速分布。李雅琪等[13]和丁雪[14]采用Fluent模拟含刚性植被的明渠流动特性。

由上述分析可知，对于含植被渠道内水流流动特性的研究大多是基于物理试验和数值模拟，并且研究流体对象都是地表上的水流，而忽视了地表下的水流。本文提出一个含植被边坡的流场计算模型，模型主要分为3层：地表上有水层和草层；地表下为土层。联立N-S方程和Biot空隙理论，推导出含刚性植被的垂向流速分布解析解。

1 理论分析

图1描绘了二维表面流在恒定情况下经过植被区的场景。在图中，很明显有3层区域：水层、草层和土壤层。

图1 含植被边坡示意Fig.1 Sketch of the vegetated slope

液体是假设为黏性的、不可压缩的均匀流，草和土层假定为均质的多孔介质，而在孔隙中运动的水视为多孔介质流。所有的物理变量和参数在一个代表截面上取平均基本体积值(REV)，从而在多孔介质内流动的流体可以用偏微分方程来描述[15]。此外，地表流和地下流假设流动已达到稳定状态。当土层达到饱和状态时，流体的垂直速度可以忽略不计。由于流动长度比流动深度大得多，意味着切向导数比法向导数小，所以切向导数可以被忽略。

1.1 基本控制方程

基于上述假设，如果水层较浅，大部分冲击动量将转移到下垫层土壤中[16]。如果忽略雨滴向水层的动量传递，同时忽略降雨在水流表层产生的水波，则Navier-Stokes方程可以简化为

x方向：

-∞

(1)

y方向：

(2)

式中：u1是水层流体在x方向(切向)上的速度，m/s；p1是水层(即层1)的流体压力，Pa；ρ和μ分别为流体的密度和动力黏滞系数，kg/m3，kg/(m·s)；g是重力加速度，m/s2；θ为坡脚，°；h1(x)为水层厚度，m；h2是草的高度，m；h(x)为水深，m；x、y为图1所示坐标。

Hsieh和Bolton[17]简化了Biot的孔隙弹性理论[18]，以模拟植被覆盖的倾斜地面上的低雷诺数流动，REV尺度多孔介质内部孔隙水的动量方程可以表示为

x方向：

(3)

y方向：

(4)

式中：uj为孔隙水在j层x方向(切向)上的速度，m/s；pj为j层孔隙介质内部孔隙水压力，Pa；μ为水的动力黏滞系数，kg/(m·s)；nj为j层孔隙度；Kpj为j层的本征渗透率，m2；下标j=2表示草层(0≤y≤h2)；下标j=3表示土层(-H≤y≤0)；H为土层厚度，m。

1.2 边界条件

在无风力、无表面水波、恒定流和使用表面压力的假设下，本研究含有以下6种边界条件：

(1) 在自由液面处(y=h(x))。

流体应力在法向(y方向)上保持连续：

p1=0

(5)

流体应力在切向(x方向)上保持连续：

(6)

(2) 在水层和草层的交界面处(y=h2)。

流速在切向(x方向)上保持连续：

u1=n2u2

(7)

流体应力在切向(x方向)上保持连续：

(8)

流体应力在法向(y方向)上保持连续：

p1=p2

(9)

(3) 在草层和土壤层的交界面处(y=0)。

流速在切向(x方向)上保持连续：

n2u2=n3u3

(10)

流体应力在切向(x方向)上保持连续：

(12)

流体应力在法向(y方向)上保持连续：

p2=p3

(13)

(4) 在土壤层底部(y=-H)。

由于无滑移条件，切向(x方向)上的流速为零：

u3=0

(13)

(5) 在边界层上部(x=0)。

u1=u2=u3=0

(14)

(6) 在边界层下部(x=L，L是沿坡地测量的坡度长度)。

边界层下部水头零深度梯度条件：

dh/dx=0

(15)

1.3 解析解

通过求解式(2)，(4)和(6)，然后将结果代入边界条件式(7)，(11)和(14)，压力分布解析解如下：

P1=ρwgcosθ(h1+h2-y)，h2

(16)

P2=ρwgcosθ(h1+h2-y)，0

(17)

P3=ρwgcosθ(h1+h2-y)，-H

(18)

式中：h1和h2分别是水层和草层的高度；θ为坡脚。

然后将压力分布代入式(1)，(3)和(5)，可得到流速分布解析解表达式：

(19)

(20)

(21)

式中：C1～C6是待求系数。

将式(19)～(21)代入式(8)～(13)和(15)，可求出C1～C6的表达式(C=n2/kp2，D=n3/kp3)。

(22)

(23)

(24)

(25)

(26)

(27)

2 分析与讨论

下文中取水的黏度μ=0.001 N·s/m2，水的密度ρw=1 000 kg/m3，重力加速度g=9.81 m/s2，坡长L=100 m。粗砂的孔隙度和渗透系数分别为n3=0.35，kp3=6.12 × 10-10m2。草层密度是由草层孔隙度决定的。植物直径(dc)与孔隙度(n2)之间的关系式及本征渗透率(kp)与孔隙率的关系式[19]为

(28)

(29)

为了验证该解的有效性，笔者发现推导的解析解在无草地面上的流速分布可以退化为Hsieh等[20]提出的解的形式(见图2)。换句话说，当坡面无植被时，Hsieh等[20]的结果是本文的一个特例。图2中的其他曲线描述了不同草密度下的速度分布。当草密度增加时，草对流体的阻力会增加，水位上升，流速降低，但草层和土壤层流速分布变得更加均匀。

图2 不同孔隙度下淹没草模型速度分布(S=0.000 1)Fig.2 Velocity distribution for submerged grass with different porosity(S=0.000 1)

草的基本参数如下：h2=0.05 m，n2=0.9，dc=0.001 m，土层厚度H=3 m。由于自然界中草的直径大小不一，且不同的直径会对坡面流体产生不同的阻力，所以图3中分别选取dc为0.003，0.005 m和0.008 m进行讨论。剪切应力可由τ=μdu/dy得到。如图3所示，其中uM和τM分别为dc=0.008 m 时的最大速度和最大剪应力。从图3(a)中，可以发现垂直速度分布中存在拐点，这与一般认知一致。从图3(b)中可知，最大剪应力存在于水草层界面或地面上，其取决于草的分布密度。此外，参考Iberall[10]，流体通过柱状纤维多孔介质时单位长度的阻力可表示为

(30)

(31)

图3 草直径对模型结果的影响Fig.3 Effects of grass diameter on the model results

图4 坡度变化对模型结果的影响Fig.4 Effects of slope ratio on the model results

土层厚度对模型流速分布的影响见图5，可知水土交界面附近有轻微的差异，总体上土层厚度对模型整体流速分布影响较小。

图5 土体厚度对模型速度分布的影响(S=0.000 1)Fig.5 Effects of the thickness of soil layer on the velocity distribution(S=0.000 1)

草的直径对阻力的影响如图6所示。流经草层的水流所受到的阻力随草径的增大而减小，而随草层孔隙度的增大而增大。这是因为相同情况下，草的直径越大孔隙就会越小。平均速度会随着孔隙度的减少而降低，从而阻力会降低。

图6 草的直径变化对阻力的影响Fig.6 Effects of different diameters of grass on drag force for various slopes

3 结 论

(1) 本文研究了水层、草层和土壤层的切向流速分布。采用N-S方程描述流体在水层的运动，通过Biot空隙理论描述流体在草层和土壤层中的运动。在速度和应力连续条件下推导出流速切向分布解析解。

(2) 通过与前人无草条件下的流速分布进行对比，验证了模型的有效性。研究发现，土层厚度变化对模型流速分布影响不大，流速分布主要受坡度的影响。

(3) 最大剪应力存在于水-草层交界面或草-土层界面上。随着草径的增大，最大剪应力主要出现在水-草层交界面处，这将有利于防治土壤侵蚀。