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

肖 毅，杜 夢，邵學軍

（清華大學水利水電工程系，水沙科學與水利水電工程國家重點實驗室，北京 100084）

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

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

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

1 二維水沙模型的改進

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

1.1 考慮邊岸植被影響的水流動量方程

根據(jù)文獻[13]推導可知正交貼體坐標下考慮植被影響的淺水動量控制方程為

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

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

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

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

2 崩岸模擬技術(shù)

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

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

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

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

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

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

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

式中：τ為平坡臨界起動應(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))，此時網(wǎng)格單元(i, j+1)參與一般水域的水流泥沙計算中，以后根據(jù)水位的變化重新確定網(wǎng)格的狀態(tài)，同時對以往的臨界陸地網(wǎng)格參數(shù)進行相應(yīng)調(diào)整．

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

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

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

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

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

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

3.1 概化游蕩型小河形成模擬試驗

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

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

表1 模擬試驗條件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ù)值模擬試驗

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

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

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

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

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

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

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

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

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

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

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

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

圖6 試驗A4小河不同時刻流場分布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)對泥沙輸運影響以及邊岸植被影響的作用，改進原有模型，并提出簡單易行的非黏性土崩岸模擬方法，成功模擬邊灘形成及演變．

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

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

［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)流橫向輸沙與河岸變形的平面二維擴展數(shù)學模型[J]. 水利學報，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］ 周 剛，王 虹，邵學軍，等. 河型轉(zhuǎn)化機理及其數(shù)值模擬(Ⅰ)：模型建立[J]. 水科學進展，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ì)河岸河道擺動的三維數(shù)值模擬[D].北京：清華大學水利系，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］ 夏軍強，王光謙.考慮河岸沖刷的彎曲河道水流及河床變形的數(shù)值模擬[J]. 水利學報，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］ 錢 寧，張 仁，周志德. 河床演變學[M]. 北京：科學出版社，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.