基于有限體積法的不同植被莖粗下坡面流水力特性數(shù)值模擬

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

（山東科技大學(xué)地球科學(xué)與工程學(xué)院，山東青島266590）

0 引言

近年來，有學(xué)者對前人的研究進行分析總結(jié)，最終將坡面流定義為以降雨為主導(dǎo)，重力為輔助動力，流經(jīng)不同坡面產(chǎn)生不同運動狀態(tài)的淺層水流［1，2］。對于現(xiàn)實自然環(huán)境而言，水力侵蝕是土壤坡面侵蝕中最為廣泛的自然現(xiàn)象（除風(fēng)力侵蝕以外）。坡面水力侵蝕的演化與水流強度、坡面土壤的抗剪切能力兩者的非協(xié)調(diào)發(fā)展機制有關(guān)，而坡面植被的存在，使得坡面水流形態(tài)更加特殊與復(fù)雜化，一方面使坡面水流內(nèi)部結(jié)構(gòu)發(fā)生改變，同時植被根系具有固結(jié)土壤作用，增強了土壤抗蝕性［3，4］，從而成為根治坡面水土流失的重要生物措施。

在國內(nèi)外學(xué)者的研究工作中，坡面流植被區(qū)域的流動特性和湍流特性受到了廣泛的關(guān)注。大量文獻提及植被分布模式對坡面流阻力的影響。Wang等［5］重點研究了植被、水流和泥沙輸移間的相互作用。植被密度和布局對流速有不同程度的影響，流阻隨著植物密度的增加而增大。Zhao 等［6］認為坡面植株的根莖能有效的降低泥沙的輸送能力，減少水土流失。張冠華等［7］通過試驗研究植被的分布格局對坡面水流阻力的影響，對不同格局的阻力系數(shù)進行對比尋找規(guī)律。Zhang 等［8］為研究植被對水流的影響，通過理論分析和實驗驗證，得出Manning糙率系數(shù)n隨雷諾數(shù)Re的增加而增大，且當植被覆蓋度越大時，Manning 糙率系數(shù)n隨雷諾數(shù)的變化率越大。姬昌輝等［9］采用水槽試驗研究了在不同水流及種植間距下的淹沒植被曼寧糙率系數(shù)的變化特征，結(jié)果表明糙率系數(shù)與植被間距呈負相關(guān)性。Ding 等［10］通過試驗得出植被分布在試驗段的下部比分布在試驗段中部及上部的攔沙效果更好，證明了植被的分布模式對明渠水流的水力特征有著重要的影響。但大都側(cè)重于實驗?zāi)M，不能準確得到植被周圍的水流流動結(jié)構(gòu)。

近年來也有一些研究人員對不同排布方式的植被進行了數(shù)值模擬的研究。Zhao和Huai［11］利用大渦模擬（LES）模型研究了明渠中不連續(xù)和淹沒植被斑塊對水流湍流的影響。Lima等［12］和Weber等［13］采用計算流體動力學(xué)（CFD）技術(shù)研究了明渠水生植被周圍的水流流動結(jié)構(gòu)。Ghani 等［14］采用數(shù)值模擬的方法，研究了圓形和交錯植被斑塊周圍的流動情況，利用不同植被密度和流速的實驗數(shù)據(jù)對數(shù)值模型進行了驗證。Anjum等［15］運用CFD 工具FLUENT，對一個縱向不連續(xù)和垂直兩層植被占據(jù)河道半寬度的水流進行了數(shù)值模擬。結(jié)果表明林隙內(nèi)的流速明顯慢于植被斑塊內(nèi)的流速，局部和不連續(xù)植被對水流結(jié)構(gòu)和阻力有較大的影響。占據(jù)半寬度的垂直分層和不連續(xù)植被的存在顯著地影響了縱向和橫向的水流。Anjum和Tanaka［16］采用k-epsilon 模型，通過覆蓋整個河道寬度的間斷垂直分層植被來闡明河道水流的流動結(jié)構(gòu)。然而這些研究大多是在相同粗度植被條件下的水流，但實際坡面植被復(fù)雜多樣，莖稈粗細程度不同，即使同一坡面上的植被，在不同的生長階段其莖稈粗度是不斷變化的，而且現(xiàn)實中洪水水流和灌溉水流在植株間的流動受植株直徑影響極大。

因此，對坡面植被植株粗度進行系統(tǒng)性的研究具有重要意義，有助于通過模型尺度對實際含植被坡面的水流阻力進行相關(guān)的預(yù)測，從而達到提高實際坡面的阻水能力，降低坡面水流的沖刷能力以及減小水力侵蝕的目的，為坡面水流阻力的研究提供價值性參考。為分析坡面植被莖桿粗度的變化對坡面流水力特性的影響，通過試驗分析研究不同流速情況下的水流阻力及流態(tài)變化規(guī)律，從而得出植被的莖桿粗度的變化對坡面流水力特性影響的定性及定量研究論點。本文的研究目的是建立一個數(shù)值模型，數(shù)值模擬是由CFD 工具FLUENT 通過采用Standardk-ε紊流模型進行的，通過均勻分布的非淹沒植被，分析坡面流的平均速度分布和湍流特性。研究還重點分析了植被間的水流分布特點。

對圓柱狀模擬植被在坡面徑流沖刷試驗和數(shù)值模擬相結(jié)合的基礎(chǔ)上，旨在研究和分析：

（1）數(shù)值模擬方法對于含植被坡面徑流沖刷試驗是否具有可行性；

（2）不同粗度變化的植被對流速的影響；

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

本數(shù)值模擬試驗以統(tǒng)一的植被行列走向與水流方向所成夾角θ= 90°、相鄰植株的行列間距a×a=60 mm×60 mm、高度為0.15 m的圓柱狀塑料棒均勻插入相應(yīng)鉆孔口徑的試驗底板。為更加直觀的反映植被莖桿粗度對坡面水動力特性的影響，本數(shù)值模擬試驗在坡度s為0 的條件下，設(shè)置3 組莖桿粗度d（3、4、5 mm）的植被，16 個單位流量q進行了模擬計算，基于有限體積法，采用Standardk-ε紊流模型，建立數(shù)值分析模型，模擬剛性植被粗度變化的條件下坡面流水流阻力的影響，對三種粗度的剛性淹沒植被明渠水流特性進行了研究。

1 數(shù)學(xué)模型

1.1 控制方程

在重力驅(qū)動下，有自由液面的坡面流動中，連續(xù)性方程為：

動量守恒方程：

式中：ui表示流速u在xi方向上的分量；ρ為流體密度；p為平均壓強；ν為運動黏度。

1.2 湍流模型

采用Standardk-ε模型［17］，不考慮流體壓縮時，Standardk-ε湍流模型的湍動能k和湍流耗散率ε的約束方程如下：

k方程：

ε方程：

式中：xj為坐標分量；ui、uj為平均相對速度分量；Pk是湍動能生成項，定義為為湍流耗散率，且ε=為動力黏滯系數(shù)；渦旋黏性系數(shù)Cμ= 0.0845。σk和σε分別為k和ε的普朗特數(shù)，σk= 1.0，σε=1.3；?；?shù)C1ε= 1.44，C2ε= 1.92，C3ε= 1.0。

2 三維湍流數(shù)值模擬分析

2.1 物理模型

本數(shù)值模擬試驗?zāi)Ｐ蛥⒄諏嶒炇以囼瀮x器建模，實體試驗儀器為長5 m，寬0.4 m，高0.3 m 的循環(huán)矩形水槽，水槽共分三個區(qū)域：上游穩(wěn)水區(qū)、試驗鋪設(shè)區(qū)、下游量水區(qū)，其中試驗鋪設(shè)區(qū)長度為3 m，上游穩(wěn)水區(qū)及下游量水區(qū)各為1 m，水槽共分為五個斷面，每個斷面之間的間距均為0.75 m。在試驗鋪設(shè)區(qū)鋪設(shè)有機樹脂板，在板中鉆孔插嵌圓柱形塑料棒模擬植被區(qū)。數(shù)值模擬試驗可以通過初始化設(shè)置達到計算所需要的水流參數(shù)條件，因此不需要單獨劃分上游穩(wěn)水區(qū)和下游量水區(qū)，所建物理模型只包含實驗鋪設(shè)區(qū)。

2.2 網(wǎng)格劃分及網(wǎng)格無關(guān)性驗證

利用ANSYS FLUENT-Mesh 繪制計算網(wǎng)格，網(wǎng)格劃分情況如圖1所示。網(wǎng)格數(shù)量達到200 多萬個，正交質(zhì)量平均值達到0.9 以上，單元質(zhì)量平均值達到了0.7 以上，表明網(wǎng)格質(zhì)量較好，滿足數(shù)值模擬運算要求。對數(shù)值計算的網(wǎng)格進行了無關(guān)性驗證（以4 mm 徑粗，0.08 m 水深為例），結(jié)果如圖2所示。隨著網(wǎng)格數(shù)量的增加，水流平均流速不斷增加，當網(wǎng)格數(shù)達到2 365 896 時，進一步增加網(wǎng)格數(shù)量，水流平均流速的增量小于0.5%，可認為此時的網(wǎng)格數(shù)已不顯著影響數(shù)值計算的結(jié)果。綜合考慮計算精度與計算效率，選用全流道網(wǎng)格總數(shù)為2 365 896的網(wǎng)格進行計算。對數(shù)值模型劃分不同數(shù)量的網(wǎng)格數(shù)并通過計算所得水流平均流速的變化來檢查網(wǎng)格數(shù)的無關(guān)性。當網(wǎng)格數(shù)對水流平均流速的變化影響小于1%時，認為網(wǎng)格數(shù)達到無關(guān)性要求。

圖1 網(wǎng)格劃分圖Fig.1 Meshing diagram

圖2 不同網(wǎng)格數(shù)下的水流平均流速Fig.2 Average flow velocity at different grid numbers

2.3 邊界條件設(shè)置

①入口邊界條件設(shè)置為質(zhì)量入口，湍流強度選擇5%，為一般水流的湍流強度；②出口邊界為壓力出口，出口壓力設(shè)置為0，溫度為300 k；③上表面利用滑移邊壁來代替水氣交界面，植被表面及壁面均采用標準壁面函數(shù)（這是更好地預(yù)測邊界層中流動的必要條件），為無滑移邊界，壁面邊界條件設(shè)置為固壁邊界。仿真和后處理在計算流體動力學(xué)（CFD）工具FLUENT 中進行。計算區(qū)域如圖1所示。采用SIMPLE 算法來實現(xiàn)壓力—速度耦合。在運算過程中，當所有殘差都低于1×10-5時，認為計算結(jié)果是收斂的。在進一步的迭代中，殘差沒有變化，并且還檢查了具有主流速的解，在后續(xù)的迭代中，該解不再明顯變化。因此，使用這些標準，可以假定解決方案已達到穩(wěn)定狀態(tài)。

3 數(shù)值模擬結(jié)果分析

3.1 3 種植被莖粗在不同流量條件下粗度對流速的影響

該數(shù)值模擬依據(jù)圖1所示模型及尺寸進行物理建模，試驗數(shù)據(jù)來自Zhang［18］等人的試驗資料，根據(jù)定床阻力試驗要求，對坡面流水力學(xué)特性進行計算機數(shù)值模擬。該模擬研究了與之前試驗數(shù)據(jù)相同流量條件下，不同植被莖粗情況下坡面流的情況下的流動狀況。模擬結(jié)果通過FLUENT后處理軟件對其進行整理分析，并與試驗數(shù)據(jù)進行對比。不同植被莖桿粗度下試驗數(shù)據(jù)與數(shù)值模擬數(shù)據(jù)水深h和平均流速v（h～v）關(guān)系曲線圖如圖3示。

圖3 3種不同植被莖稈粗度h～v關(guān)系曲線圖Fig.3 Three different vegetation stalk thickness h～v curve diagram

根據(jù)數(shù)值模擬結(jié)果發(fā)現(xiàn)，3種不同植被莖稈粗度的h～v關(guān)系總體呈現(xiàn)出隨著流量的不斷增大，水深隨之增加，流速不斷增大；當水深h較小時，流速隨水深的變化較大，且不同植被莖稈粗度的速度差異較大，隨著水深的增加，流速呈緩慢變化趨勢。在同一水深h下，不同植被莖稈粗度情況下，水流流速不同；植被莖稈粗度越小對應(yīng)的流速值也越大，3 種植被莖稈粗度的水流流速關(guān)系為：v3mm＞v4mm＞v5mm，這一規(guī)律說明坡面植被莖稈粗度的微小變化對水流流速有一定的影響，植被莖稈粗度越大，水流流速越小。

試驗結(jié)果與數(shù)值模擬結(jié)果進行對比如圖3 示，數(shù)值模擬結(jié)果比試驗實測數(shù)據(jù)偏大，但保持在10%以內(nèi)，因此可以認為數(shù)值模擬在坡面流水流阻力試驗研究中的應(yīng)用是可行的。對于坡面流而言，影響阻力系數(shù)的因素眾多，在本次對比驗證試驗中，造成這種差異的原因可能是：在數(shù)值模擬試驗中，雨滴打擊力在試驗中未涉及，因而可忽略不計；數(shù)值模擬試驗中重力和水壓力也已通過參數(shù)設(shè)置完成，水槽壁面條件設(shè)置為光滑壁面，為理想環(huán)境。從以上分析可見，誤差產(chǎn)生的主要原因可確定為坡面流表面張力的作用以及光滑的壁面條件。

3.2 3 種植被莖粗在相同流量條件下水流時均流速分布

本數(shù)值模擬選取3、4、5 mm 3 種植被莖粗分為3 組，計算了在相同流量和水深植被淹沒狀態(tài)下的水流流速。并且每組模擬中各選取了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 組數(shù)據(jù)，總共采集了2 400 個數(shù)據(jù)，保證了垂向速度分布的準確性。

圖4 計算域示意圖（單位：mm）Fig.4 Schematic diagram of calculation domain

3.2.1 縱向流速分布

圖5～7所示分別為3、4、5 mm 直徑的植被在相同流量條件下各測點水流流速的垂向分布?？傮w看來，處在兩行植被中間的測點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 位于植被正下游且緊鄰植被，受到植被阻力影響最大，其縱向流速變化最小，相對于其他測點的縱向流速小的多。這是因為水流經(jīng)過植被時做圓柱繞流，水流流經(jīng)植被前受到植被的阻擋，迎流面壓強勢能逐漸增大，動能逐漸變小，因而速度逐漸變小。測點a 縱向流速次之，測點c 和e 速度相近，這是因為測點c 和e 均位于兩植被之間，離植被較遠，水流受到植被阻力影響較小，且由于兩測點相距較近，可以忽略流速損失。處在兩行植被中間的測點l～p中，測點n縱向流速最大，達到0.298 4 m/s，且為所有測點流速中最大。通過（圖8）測點i、nyx剖面流速分布云圖也可以看到這一現(xiàn)象，測點n 位于兩株植被中間，流速表現(xiàn)最大。測點i所處位置更接近植被，流速表現(xiàn)小于測點n。我們還可以看出，從兩株植被中間位置（測點n）到植被位置，流速表現(xiàn)為逐漸變小。這與植被的擋水作用有關(guā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 垂向流速分布

考慮到植被對水流流動結(jié)構(gòu)的影響較大，有必要對含植被水流的垂向流速變化進行研究。從3種粗度植被測點縱向流速的垂向分布（圖5～7）可以很直觀的看出，在模擬的植被區(qū)域內(nèi)，速度的垂直變化是單調(diào)的，最小速度發(fā)生在河床。這與Zeng C［19］的測量結(jié)果相同。與Neumeier［20］的測量結(jié)果不同，Neumeier 的測量結(jié)果表明，最小流速發(fā)生在地層上方約5 cm處。主要原因是Neumeir實驗中采用了變截面自然植被。最大植被密度出現(xiàn)在床上5 cm 左右，對應(yīng)最小流速水平。位于不同xz剖面的測點a～e、f～k、l～p垂向流速分布呈現(xiàn)很明顯的規(guī)律，垂向流速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剖面垂向流速分布云圖可以直觀地體現(xiàn)出流速變化。

以（圖5）3 mm 直徑植被測點縱向流速的垂向分布為例，測點f～k 垂向流速分布幾乎一致，植被淹沒高度h=0.02 m 左右時出現(xiàn)拐點，垂向流速趨近于最大值（0.219 8 ～0.239 9 m/s）。測點l～p 的垂向流速分布一致性更強，植被淹沒高度h=0.03 m 左右時出現(xiàn)拐點，縱向流速趨近于最大值（0.287 3 ～0.298 4 m/s），它們的流速垂向分布曲線近乎重合。這一現(xiàn)象通過（圖10）測點a-pxz剖面垂向流速分布云圖能夠直觀體現(xiàn)。由圖10（b）和（c）可以發(fā)現(xiàn)f～k、l～p 測點所在xz剖面的垂向速度分布較為均勻，所以各個測點的流速垂向分布相差較小。與之相比，測點a～e 所在xz剖面的垂向速度分布變化較大，植被淹沒高度h=0.02 m 左右時出現(xiàn)拐點，垂向流速趨近于最大值（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 垂向流速變化幅度最?。ㄆ渲袦y點d 變化幅度最大，為0～0.147 9 m/s）?？拷脖坏臏y點f～k 垂向流速變化幅度大于測點a～e（其中測點i 變化幅度最大，為0 ～0.239 9 m/s）。處在兩行植被中間的測點l～p 的垂向流速變化幅度最大（其中測點n 變化幅度最大，為0～0.298 4 m/s）。通過（圖9）測點a～pyz剖面垂向流速分布云圖可以看出這種規(guī)律。且3組粗度植被的測點縱向流速的垂向分布圖都呈現(xiàn)這種規(guī)律。通過進一步分析發(fā)現(xiàn)，這與植被對水流的紊動作用有關(guān)［21］。測點a～e 與植被同行分布，由于植被與水流交界面處水流紊動較為劇烈，紊動強度明顯增加，其紊流強度垂向分布，在植物頂部的水流劇烈運動產(chǎn)生較大的能量損失是形成水流阻力的重要原因，因此測點a～e 垂向流速變化幅度最小。而測點f～k 與植被排布處較遠，植被對水流的紊動作用減小，因此測點f～k 垂向流速變化幅度大于測點a～e。處在兩行植被中間的測點l～p 離植被分布區(qū)最遠，植被對水流的紊動作用最小，因而垂向流速變化幅度最大。

圖9 測點（a～p）yz剖面垂向流速分布云圖（單位：m/s）Fig.9 Measurement point（a～p）yz profile vertical velocity distribution cloud map

另外需要注意的是，在數(shù)值模擬數(shù)據(jù)結(jié)果中，在水槽底部有速度突增，Liu等［22］，Huai等［23］和Pang等［24］在植被水流的試驗中也發(fā)現(xiàn)了同樣的現(xiàn)象，Liu 等［22］認為水槽底部的速度突增可能是由植被底部的馬蹄形渦旋或是結(jié)合渦旋引起的，馬蹄形渦旋會使得水流從周圍區(qū)域較快的繞過植被底部，從而導(dǎo)致植被底部附近流速突增。

4 結(jié)論

采用標準紊流模型，建立數(shù)值分析模型，進行三維數(shù)值計算，通過在流場入口設(shè)置為質(zhì)量入口的方式，對矩形明渠中不同莖徑連續(xù)植被水流的流動進行了模擬。計算所得的流速值與實測值吻合良好，驗證了數(shù)值模擬方法的正確性和可行性。在建立數(shù)值分析模型的基礎(chǔ)上，對坡面水流的流動過程中的水動力學(xué)特性的變化規(guī)律進行分析，可以得出以下結(jié)論。

（1）試驗結(jié)果與數(shù)值模擬結(jié)果進行對比分析發(fā)現(xiàn)，模擬結(jié)果比試驗實測數(shù)據(jù)偏大，但誤差保持在10%以內(nèi)，與試驗值吻合良好，驗證了數(shù)值模擬在坡面流水流阻力試驗研究中的可行性。

（2）根據(jù)數(shù)值模擬結(jié)果發(fā)現(xiàn)，3 種不同植被莖稈粗度的h～v關(guān)系總體呈現(xiàn)出隨著流量的不斷增大，水深隨之增加，流速不斷增大；植被莖桿粗度增大，流速減小。這是由于植被莖桿粗度增大，相對應(yīng)的總體覆蓋度增大，水流與植被接觸面增大，從而使水流阻力增大

（3）縱向流速分布規(guī)律：總體看來，處在兩行植被中間的測點l～p 縱向流速相對較大，靠近植被的測點f-k 縱向流速次之，與植被同行分布的測點a～e縱向流速最小。水流縱向流速隨著到植被距離的增加先增加后減小。且隨著植被粗度的增加，水流縱向流速逐漸變小。這表明植被根部附近有利于沉積物的積聚，對植物群落以及水生生物提供營養(yǎng)物質(zhì)環(huán)境和空氣。

（4）垂向流速分布規(guī)律：在y方向距植被越遠，垂向流速的增長率逐漸減小。從3 種粗度植被垂向分布圖5～7 可以很直觀的看出，位于不同xz剖面的測點a～e、f～k、l～p 垂向流速分布呈現(xiàn)很明顯的規(guī)律，垂向流速v（測點l～p）＞垂向流速v（測點f～k）＞垂向流速v（測點a～e），且流速分布差異非常明顯。測點a～e垂向流速變化幅度最小，測點f-k 垂向流速變化幅度大于測點a～e。處在兩行植被中間的測點l～p 離植被分布區(qū)最遠，植被對水流的紊動作用最小，因而垂向流速變化幅度最大。植被間隙的低流速有利于泥沙沉積物沉積，有助于坡面穩(wěn)定，維持生態(tài)環(huán)境?！?/p>