水-草-土耦合邊坡流場特性研究

陳 彬，姚 錫 偉，毛 朝 勛

(內江職業技術學院，四川 內江 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) 最大剪應力存在于水-草層交界面或草-土層界面上。隨著草徑的增大，最大剪應力主要出現在水-草層交界面處，這將有利于防治土壤侵蝕。