肖 毅，杜 夢(mèng)，邵學(xué)軍

（清華大學(xué)水利水電工程系，水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084）

游蕩河型是天然河流中常見的河型之一，在國(guó)內(nèi)外分布較廣．我國(guó)黃河中下游均為典型的游蕩型河道，因此采用多種手段研究游蕩型河流的演變規(guī)律，在理論與實(shí)踐上都有重要意義．20世紀(jì)80年代以來(lái)，隨著計(jì)算數(shù)學(xué)的不斷發(fā)展，數(shù)值模擬已作為重要手段應(yīng)用于河型成因機(jī)理研究中．基于其簡(jiǎn)單性、易操作性等優(yōu)點(diǎn)，文獻(xiàn)[1-8]分別通過(guò)改進(jìn)水沙動(dòng)力學(xué)模型，模擬彎曲河道形成及發(fā)展過(guò)程，研究其形成機(jī)理及各因素之間的關(guān)系．

相較于彎曲河道模擬研究的廣泛性，目前對(duì)游蕩河道的數(shù)值研究及討論仍然較為缺乏，Murray等[9]通過(guò)自動(dòng)元胞機(jī)得到辮狀河道的流場(chǎng)；夏軍強(qiáng)等[10]通過(guò)改進(jìn)二維崩岸模塊模擬黃河下游游蕩河段的發(fā)展；Schuurman等[11]通過(guò)對(duì)比二維及三維數(shù)值模型對(duì)游蕩小河的模擬，研究其沙洲形成及演變過(guò)程，文獻(xiàn)[12]將室內(nèi)辮狀小河形成過(guò)程進(jìn)行模擬，得到結(jié)果與試驗(yàn)數(shù)據(jù)基本吻合．但上述河型數(shù)值模擬的研究均局限于對(duì)單一河型形成過(guò)程的模擬，缺乏對(duì)某一河型形成及轉(zhuǎn)化過(guò)程的完整的模擬．

筆者基于改進(jìn)的平面二維數(shù)值模型，建立概化河道的數(shù)值試驗(yàn)，模擬游蕩型小河的形成；在此基礎(chǔ)上，通過(guò)改變來(lái)水量、來(lái)沙量及邊岸植被等控制因素，分別模擬游蕩小河進(jìn)一步發(fā)展；分析討論上述因素對(duì)游蕩河型轉(zhuǎn)化的影響，并將結(jié)論與現(xiàn)有河型轉(zhuǎn)化理論對(duì)比，說(shuō)明利用數(shù)值模擬技術(shù)研究游蕩河型轉(zhuǎn)化及成因的可行性．

1 二維水沙模型的改進(jìn)

本文采用的基本平面二維水沙模型詳見文獻(xiàn)[7]．

1.1 考慮邊岸植被影響的水流動(dòng)量方程

根據(jù)文獻(xiàn)[13]推導(dǎo)可知正交貼體坐標(biāo)下考慮植被影響的淺水動(dòng)量控制方程為

1.2 考慮河床結(jié)構(gòu)及河灣形態(tài)對(duì)泥沙輸運(yùn)的影響

考慮河床結(jié)構(gòu)對(duì)泥沙輸運(yùn)影響，本節(jié)采用泥沙輸運(yùn)方向[14]為

式中：n為河床糙率；vavg為該點(diǎn)平均流速；γ 為水的容重；γs為泥沙的容重．

由于現(xiàn)有二維平面模型無(wú)法準(zhǔn)確模擬彎道二次流對(duì)泥沙輸運(yùn)作用，為此對(duì)近底切應(yīng)力方向進(jìn)行修正，其表達(dá)式為

2 崩岸模擬技術(shù)

本文基于非黏性土力學(xué)崩岸機(jī)理，在不劃分新網(wǎng)格的基礎(chǔ)上，只增加崩岸記憶數(shù)組，提高計(jì)算效率．可以成功模擬邊岸崩塌及沙洲形成過(guò)程，具體模擬方法如下．

2.1 考慮河岸形態(tài)對(duì)崩岸的影響

影響坍岸強(qiáng)度的主要因素是水流強(qiáng)度，河岸土質(zhì)條件和河灣形態(tài)[15]．tΔ時(shí)間內(nèi)橫向沖刷距離根據(jù)經(jīng)驗(yàn)公式[16]得

式中：WΔ為河岸tΔ時(shí)間內(nèi)的橫向后退距離：lC為河岸橫向沖刷系數(shù)，主要取決于河岸土體的物理化學(xué)特性，根據(jù)實(shí)際情況進(jìn)行率定；τ為邊岸水流切應(yīng)力；cτ為邊岸土體發(fā)生沖刷的起動(dòng)切應(yīng)力；bkγ為邊岸土體的容重．

研究表明[17]，凹岸沖刷的強(qiáng)弱與縱向切應(yīng)力成正比，而在彎道水流中，床面切應(yīng)力的分布與縱向流速分布一致．基于原模型對(duì)彎道二次流影響的考慮，修正近底切應(yīng)力可得

式中：au為斷面平均流速；cu為縱向垂線平均流速．

河岸泥沙與河床泥沙的受力特性基本一致，但是由于河岸邊灘具有一定坡度，甚至處于直立狀態(tài)，導(dǎo)致其起動(dòng)情況發(fā)生改變，因此需要修正泥沙的臨界起動(dòng)切應(yīng)力，其表達(dá)式[1]為

式中：τ為平坡臨界起動(dòng)應(yīng)力，采用 Leo C van Rijn公式；β為縱向坡腳；α為河岸岸坡；φ為泥沙水下休止角．

2.2 河岸穩(wěn)定性判定及數(shù)值模擬

本模型采用崩岸數(shù)值模擬方法如圖1所示．

由圖1有

ΔB+ΔL＞ x( i, j + 1 ) - x (i, j )(圖 1(a))，此時(shí)網(wǎng)格單元(i, j+1)參與一般水域的水流泥沙計(jì)算中，以后根據(jù)水位的變化重新確定網(wǎng)格的狀態(tài)，同時(shí)對(duì)以往的臨界陸地網(wǎng)格參數(shù)進(jìn)行相應(yīng)調(diào)整．

ΔB+ΔL≤ x( i,j) -x ( i-1 ,j) ≤ΔB+ΔL+Δp (圖1(b))，此時(shí)網(wǎng)格單元(i, j + 1 )地形根據(jù)幾何關(guān)系成比例進(jìn)行調(diào)整，并同時(shí)采用記憶數(shù)組網(wǎng)格記錄該方向的崩岸信息，納入下一次的崩岸計(jì)算．

x( i, j+1)-x(i,j)＞ΔB+ΔL+Δp (圖 1(c))，此時(shí)網(wǎng)格單元 (i, j + 1 )不進(jìn)行調(diào)整，采用該方向的記憶單元，記錄該方向的變形，納入下一次的崩岸計(jì)算中．

本文采用計(jì)算臨界水域點(diǎn) 1相鄰網(wǎng)格最大剪切應(yīng)力法來(lái)判定其崩塌方向(圖 1(d))；不再局限于固定河岸邊界，使其能夠模擬沙洲邊灘的形成．根據(jù)橫向崩塌后沙量守恒原則，將崩岸后的土體按比例平鋪于相鄰的計(jì)算網(wǎng)格節(jié)點(diǎn)上，調(diào)整各相鄰點(diǎn)的地形．

3 游蕩型河流轉(zhuǎn)化控制因素的數(shù)值模擬

游蕩型河流演變是多因素相互作用的復(fù)雜響應(yīng)過(guò)程，由于試驗(yàn)的限制和實(shí)測(cè)資料的缺乏，對(duì)其河型轉(zhuǎn)化的控制因素的研究仍未得到系統(tǒng)的結(jié)論；本節(jié)將首先利用改進(jìn)模型模擬游蕩型小河的形成；然后通過(guò)調(diào)整其來(lái)水量、來(lái)沙量及邊岸植被控制因素，進(jìn)一步模擬該小河平面形態(tài)變化．通過(guò)模擬結(jié)果分析上述控制因素對(duì)游蕩型小河轉(zhuǎn)化的影響，并將結(jié)論與現(xiàn)有游蕩河型轉(zhuǎn)化理論進(jìn)行對(duì)比．

圖1 非黏性土崩岸計(jì)算方法Fig.1 Method for calculation of non-cohesive bank erosion

3.1 概化游蕩型小河形成模擬試驗(yàn)

概化河道長(zhǎng) 10,000,m，寬 300,m，河道比降為0.4×10-3．其初始地形高程見圖 2，模擬條件見表 1，河岸組成與床沙相同，均為非黏性均勻沙．計(jì)算過(guò)程中，水流時(shí)間步長(zhǎng)為 6,s，泥沙時(shí)間步長(zhǎng)為 12,s，模擬總歷時(shí) 760,d，實(shí)際計(jì)算用時(shí)約 20,h．模擬小河地形高程見圖3．

從圖3可知，模擬720,d后概化河道形成了游蕩散亂河型，由模擬條件表 1可得，河道在 360,d后來(lái)水來(lái)沙量同時(shí)增加，從而導(dǎo)致了水流攜帶泥沙能力的增強(qiáng)，從上游攜帶大量泥沙，隨著河道展寬，其挾沙能力降低，大量的泥沙淤積在下游河床上，導(dǎo)致水流散亂分汊，形成了游蕩型河道．

表1 模擬試驗(yàn)條件Tab.1 Simulation of experimental conditions

圖2 概化河道初始平面地形Fig.2 Layout of conceptual river channel

圖3 概化河道720,d后平面地形Fig.3 Layout of the conceptual river channel after 720,days

3.2 游蕩型小河轉(zhuǎn)化控制因素影響的數(shù)值模擬試驗(yàn)

基于上述產(chǎn)生的游蕩型小河的地形，本小節(jié)通過(guò)概化數(shù)值試驗(yàn)，改變控制因素，模擬其河道進(jìn)一步平面變形，試驗(yàn)分為 4組，試驗(yàn)條件及模擬結(jié)果見表2．計(jì)算過(guò)程中，水流時(shí)間步長(zhǎng)為 6,s，泥沙時(shí)間步長(zhǎng)為12,s，試驗(yàn)小河模擬歷時(shí)600,d后地形高程見圖4．

表2 數(shù)值試驗(yàn)條件及模擬結(jié)果Tab.2 Experimental conditions and simulation results

圖4 試驗(yàn)小河600,d后平面形態(tài)Fig.4 Layout of experimental channel after 600,days

3.3 各控制因素對(duì)游蕩河型轉(zhuǎn)化影響分析

3.3.1 單一控制因素對(duì)游蕩型小河平面形態(tài)變化的影響分析

圖 4(a)～(c)分別為試驗(yàn) A1～A3減小流量，減小來(lái)沙量及邊岸植被覆蓋的單因素改變下游蕩型小河 600,d后的平面形態(tài)．試驗(yàn) A1中減小流量相當(dāng)于小水帶大沙，此時(shí)的水流挾沙力小于其含沙量，因此河道上游的部分汊道將被淤死，出現(xiàn)了主槽；而在下游大量泥沙淤積于床面，使河道仍然散亂(圖 4(a))；試驗(yàn) A2減小來(lái)沙量，相當(dāng)于大水帶小沙，此時(shí)的河道挾沙力大于其含沙量，河道將沖深主槽，由于大水趨直，因此河道的下游出現(xiàn)了一段順直河段，且汊道減少(圖 4(b))；試驗(yàn) A3中加入邊岸植被影響，其邊岸的抗沖性增強(qiáng)，限制河道進(jìn)一步的展寬，由于邊灘植草起到穩(wěn)定邊灘的作用，因此河型形態(tài)變化不大，一些汊道由于沖深而使水流歸槽，從而其散亂程度有所降低(圖4(c))．

選取試驗(yàn) A1～A3組小河 6,000,m 處斷面形態(tài)與游蕩小河初始斷面形態(tài)進(jìn)行比較(圖 5)：試驗(yàn) A1的斷面形態(tài)仍然為多汊道，但是其汊道數(shù)較初始為少；試驗(yàn) A2的斷面形態(tài)已經(jīng)由原來(lái)的多汊道的“W”型轉(zhuǎn)變?yōu)橛兄鞑鄣摹癠”形；試驗(yàn) A3小河由于邊岸植被固岸的影響，因此展寬小于前兩組，其形態(tài)與初始形態(tài)基本一致，主槽的刷深較初始地形加深，說(shuō)明邊岸植被有維持原有河道形態(tài)的作用．

圖5 試驗(yàn)小河600 d后6,000 m斷面形態(tài)對(duì)比Fig.5 Comparison of bed deformation at the 6,000,m sections of channel after 600 days

3.3.2 控制因素綜合作用對(duì)游蕩型小河平面形態(tài)變化的影響分析

圖 4(d)為試驗(yàn) A4控制變量(來(lái)水量，來(lái)沙量，邊岸植被覆蓋)綜合作用下模擬小河 600,d后的平面形態(tài)．從表 2試驗(yàn)條件可知，隨著來(lái)水來(lái)沙的減小，游蕩型小河床面上的泥沙淤積減小，散亂程度減小，由于小水下水流歸槽的作用，使得試驗(yàn)小河逐漸形成主槽；由于此時(shí)邊岸植被固岸的共同作用下，小河的主槽能夠得到進(jìn)一步的維持，使其轉(zhuǎn)化為彎曲小河的速度加快；當(dāng)不考慮邊岸植被影響只考慮來(lái)水來(lái)沙量減少時(shí)的游蕩型小河轉(zhuǎn)化為彎曲型小河需要的時(shí)間為 500,d，說(shuō)明邊岸植被起到穩(wěn)定主槽作用，促進(jìn)其游蕩型小河向彎曲型小河轉(zhuǎn)化．圖6為試驗(yàn)A4模擬小河不同時(shí)刻的流場(chǎng)，可以看到來(lái)水來(lái)沙的減少，使水流逐漸歸于主槽形成單一的河槽，從而促使其向彎曲型小河轉(zhuǎn)化．

3.3.3 試驗(yàn)結(jié)論與現(xiàn)有游蕩型河流轉(zhuǎn)化理論對(duì)比

根據(jù)表 3[15]可知，當(dāng)流量、來(lái)沙量減小，植被覆蓋率增加的情況下，河型發(fā)展趨勢(shì)為游蕩向彎曲小河轉(zhuǎn)化，通過(guò)數(shù)值試驗(yàn) A1～A3可以得到：當(dāng)僅改變單一控制因素時(shí)，游蕩型小河的轉(zhuǎn)化速度明顯小于控制因素綜合作用下的轉(zhuǎn)化過(guò)程(試驗(yàn) A4)；而通過(guò)試驗(yàn)A3對(duì)植被覆蓋率對(duì)河型轉(zhuǎn)化的影響模擬中，當(dāng)其原河道較為散亂時(shí)，增加邊岸植被可以穩(wěn)定現(xiàn)有的河道，但是對(duì)其轉(zhuǎn)化為彎曲河道的影響并不十分明顯，這可能與模擬歷時(shí)不夠有關(guān)，但是在模擬結(jié)果的平面形態(tài)圖4(c)中仍可以看到其汊道逐漸減少．

圖6 試驗(yàn)A4小河不同時(shí)刻流場(chǎng)分布Fig.6 Temporal changes in flow field including velocity and surface elevation of test A4

表3 控制條件發(fā)生改變河型的發(fā)展方向[15]Tab.3 Transformation of channel with change of control factors[15]

4 結(jié) 論

(1) 考慮彎道二次流、河床結(jié)構(gòu)對(duì)泥沙輸運(yùn)影響以及邊岸植被影響的作用，改進(jìn)原有模型，并提出簡(jiǎn)單易行的非黏性土崩岸模擬方法，成功模擬邊灘形成及演變．

(2) 基于改進(jìn)模型，建立概化游蕩型小河轉(zhuǎn)化控制因素影響的數(shù)值試驗(yàn)，首先模擬游蕩型小河的形成；在此基礎(chǔ)上，通過(guò)調(diào)整來(lái)水來(lái)沙量及邊岸植被覆蓋率等控制因素，分別模擬游蕩小河進(jìn)一步發(fā)展；通過(guò)模擬結(jié)果對(duì)控制游蕩河型轉(zhuǎn)化的因素進(jìn)行分析得到：大水小沙、小水小沙、以及邊岸植被覆蓋率的增加都有使游蕩型小河變?yōu)閺澢托『拥内厔?shì)；小水小沙以及加強(qiáng)邊岸的植被覆蓋率可以加速?gòu)澢『拥男纬?；結(jié)論與現(xiàn)有河型轉(zhuǎn)化理論吻合，說(shuō)明利用數(shù)值模擬技術(shù)研究游蕩河型轉(zhuǎn)化及成因的可行性．

(3) 由于本模型考慮彎道二次流、河床結(jié)構(gòu)及植被影響的公式與參數(shù)都是基于理論與試驗(yàn)相結(jié)合得到經(jīng)驗(yàn)公式，因此泥沙的輸運(yùn)過(guò)程以及河道的橫向崩岸過(guò)程將受到影響，從而導(dǎo)致河道平面形態(tài)模擬的不準(zhǔn)確(如河道崩岸系數(shù)過(guò)大而導(dǎo)致河岸展寬過(guò)快造成流速降低，從而影響泥沙的輸運(yùn)及其進(jìn)一步的形態(tài)變化)；同時(shí)，由于采用數(shù)組記憶邊岸形態(tài)的變化過(guò)程，不能及時(shí)反映河岸崩塌過(guò)程對(duì)水流之間的相互關(guān)系．因此如何找到適應(yīng)范圍更廣、理論上更合理的公式并改進(jìn)數(shù)值計(jì)算方法，更為準(zhǔn)確反映天然河道的演變過(guò)程將是未來(lái)河道演變模擬研究的重點(diǎn)．

［1］ Duan J G，Julien P Y. Numerical simulation of the inception of channel meandering[J].Earth Surface Processes and Landforms，2005，30(7)：1093-1110.

［2］ Darby S E，Alabyan A M，Van de Wiel M J.Numerical simulation of bank erosion and channel migration in meandering rivers[J].Water Resources Research，2002，38(9)：1-21.

［3］ Olsen N R B. 3D CFD modeling of a self-forming meandering channel [J].Journal of Hydraulic Engineering，2003，129(5)：366-372.

［4］ Ruther N，Olsen N R B. CFD modeling of meandering river evolution[C]//Proceedings of30th IAHR Congress.Thessaloniki，Greece，2003.

［5］ Dulal K P，Kobayashi K，Shimizu Y，et al. Numerical computation of free meandering channels with the application of slump blocks on the outer bends[J].Journal of Hydro-Environment Research，2010，3(4)：239-246.

［6］ 鐘德鈺，張洪武. 考慮環(huán)流橫向輸沙與河岸變形的平面二維擴(kuò)展數(shù)學(xué)模型[J]. 水利學(xué)報(bào)，2004(7)：14-20.

Zhong Deyu，Zhang Hongwu. Extended 2-D numerical model for alluvial river considering transverse transport of sediment and bank erosion due to secondary flow in river bends[J].Journal of Hydraulic Engineering，2004(7)：14-20(in Chinese).

［7］ 周 剛，王 虹，邵學(xué)軍，等. 河型轉(zhuǎn)化機(jī)理及其數(shù)值模擬(Ⅰ)：模型建立[J]. 水科學(xué)進(jìn)展，2010(3)：145-152.

Zhou Gang，Wang Hong，Shao Xuejun，et al. Mechanism of channel pattern changes and its numerical simulation(Ⅰ)：Numerical model[J].Advances in Water Science，2010(3)：145-152(in Chinese).

［8］ 假冬冬．非均質(zhì)河岸河道擺動(dòng)的三維數(shù)值模擬[D].北京：清華大學(xué)水利系，2010.

Jia Dongdong. 3-D Numerical Simulation of Lateral Migration of Alluvial Channels with Composite Banks [D].Beijing：Department of Hydraulic Engineering，Tsinghua University，2010(in Chinese).

［9］ Murray A B，Paola C. A cellular model of braided rivers[J].Nature，1994，371(6492)：54-57.

［10］ 夏軍強(qiáng)，王光謙.考慮河岸沖刷的彎曲河道水流及河床變形的數(shù)值模擬[J]. 水利學(xué)報(bào)，2002(6)：60-66.

Xia Junqiang，Wang Guangqian. Numerical simulation of flow and bed deformation in meandering rivers considering the erosion of bank[J].Journal of Hydraulic Engineering，2002(6)：60-66(in Chinese).

［11］ Schuurman F，Kleinhans M. Self-formed braided bar pattern in a numerical model[J].River，Coastal and Estuarine Morphodynamics：RCEM2011.

［12］ Yang Haiyan，Lin Binliang，Zhou Jianjun. Braided river modelling using a numerical model[J].River，Coastal and Estuarine Morphodynamics：RCEM2011.

［13］ Ikeda S，Izumi N. Width and depth of self-formed straight gravel rivers with bank vegetation[J].Water Resources Research，1990，26(10)：2353-2364.

［14］ Koch F G，F(xiàn)lokstra C. Bed level computations for curved alluvial channels[C]//Proceedings of theXIXth Congress of the IAHR. New Delhi，India，1981：357-364.

［15］ 錢 寧，張 仁，周志德. 河床演變學(xué)[M]. 北京：科學(xué)出版社，1987.

Chien Ning，Zhang Ren，Zhou Zhide.Fluvial Processes[M]. Beijing：Science Press，1987(in Chinese).

［16］ Osman A M，Thorne C R. Riverbank stability analysis(Ⅰ)：Theory[J].Journal of Hydraulic Engineering，1988，114(2)：134-150.

［17］ Varshney D V，Garde R J. Shear distribution in bend in rectangular channels[J].Journal of the Hydraulics Division，ASCE，1975，101(8)：1053-1066.