游荡型河流转化控制因素的二维数值模拟

肖 毅，杜 梦，邵学军

（清华大学水利水电工程系，水沙科学与水利水电工程国家重点实验室，北京 100084）

游荡河型是天然河流中常见的河型之一，在国内外分布较广．我国黄河中下游均为典型的游荡型河道，因此采用多种手段研究游荡型河流的演变规律，在理论与实践上都有重要意义．20世纪80年代以来，随着计算数学的不断发展，数值模拟已作为重要手段应用于河型成因机理研究中．基于其简单性、易操作性等优点，文献[1-8]分别通过改进水沙动力学模型，模拟弯曲河道形成及发展过程，研究其形成机理及各因素之间的关系．

相较于弯曲河道模拟研究的广泛性，目前对游荡河道的数值研究及讨论仍然较为缺乏，Murray等[9]通过自动元胞机得到辫状河道的流场；夏军强等[10]通过改进二维崩岸模块模拟黄河下游游荡河段的发展；Schuurman等[11]通过对比二维及三维数值模型对游荡小河的模拟，研究其沙洲形成及演变过程，文献[12]将室内辫状小河形成过程进行模拟，得到结果与试验数据基本吻合．但上述河型数值模拟的研究均局限于对单一河型形成过程的模拟，缺乏对某一河型形成及转化过程的完整的模拟．

笔者基于改进的平面二维数值模型，建立概化河道的数值试验，模拟游荡型小河的形成；在此基础上，通过改变来水量、来沙量及边岸植被等控制因素，分别模拟游荡小河进一步发展；分析讨论上述因素对游荡河型转化的影响，并将结论与现有河型转化理论对比，说明利用数值模拟技术研究游荡河型转化及成因的可行性．

1 二维水沙模型的改进

本文采用的基本平面二维水沙模型详见文献[7]．

1.1 考虑边岸植被影响的水流动量方程

根据文献[13]推导可知正交贴体坐标下考虑植被影响的浅水动量控制方程为

1.2 考虑河床结构及河湾形态对泥沙输运的影响

考虑河床结构对泥沙输运影响，本节采用泥沙输运方向[14]为

式中：n为河床糙率；vavg为该点平均流速；γ 为水的容重；γs为泥沙的容重．

由于现有二维平面模型无法准确模拟弯道二次流对泥沙输运作用，为此对近底切应力方向进行修正，其表达式为

2 崩岸模拟技术

本文基于非黏性土力学崩岸机理，在不划分新网格的基础上，只增加崩岸记忆数组，提高计算效率．可以成功模拟边岸崩塌及沙洲形成过程，具体模拟方法如下．

2.1 考虑河岸形态对崩岸的影响

影响坍岸强度的主要因素是水流强度，河岸土质条件和河湾形态[15]．tΔ时间内横向冲刷距离根据经验公式[16]得

式中：WΔ为河岸tΔ时间内的横向后退距离：lC为河岸横向冲刷系数，主要取决于河岸土体的物理化学特性，根据实际情况进行率定；τ为边岸水流切应力；cτ为边岸土体发生冲刷的起动切应力；bkγ为边岸土体的容重．

研究表明[17]，凹岸冲刷的强弱与纵向切应力成正比，而在弯道水流中，床面切应力的分布与纵向流速分布一致．基于原模型对弯道二次流影响的考虑，修正近底切应力可得

式中：au为断面平均流速；cu为纵向垂线平均流速．

河岸泥沙与河床泥沙的受力特性基本一致，但是由于河岸边滩具有一定坡度，甚至处于直立状态，导致其起动情况发生改变，因此需要修正泥沙的临界起动切应力，其表达式[1]为

式中：τ为平坡临界起动应力，采用 Leo C van Rijn公式；β为纵向坡脚；α为河岸岸坡；φ为泥沙水下休止角．

2.2 河岸稳定性判定及数值模拟

本模型采用崩岸数值模拟方法如图1所示．

由图1有

ΔB+ΔL＞ x( i, j + 1 ) - x (i, j )(图 1(a))，此时网格单元(i, j+1)参与一般水域的水流泥沙计算中，以后根据水位的变化重新确定网格的状态，同时对以往的临界陆地网格参数进行相应调整．

ΔB+ΔL≤ x( i,j) -x ( i-1 ,j) ≤ΔB+ΔL+Δp (图1(b))，此时网格单元(i, j + 1 )地形根据几何关系成比例进行调整，并同时采用记忆数组网格记录该方向的崩岸信息，纳入下一次的崩岸计算．

x( i, j+1)-x(i,j)＞ΔB+ΔL+Δp (图 1(c))，此时网格单元 (i, j + 1 )不进行调整，采用该方向的记忆单元，记录该方向的变形，纳入下一次的崩岸计算中．

本文采用计算临界水域点 1相邻网格最大剪切应力法来判定其崩塌方向(图 1(d))；不再局限于固定河岸边界，使其能够模拟沙洲边滩的形成．根据横向崩塌后沙量守恒原则，将崩岸后的土体按比例平铺于相邻的计算网格节点上，调整各相邻点的地形．

3 游荡型河流转化控制因素的数值模拟

游荡型河流演变是多因素相互作用的复杂响应过程，由于试验的限制和实测资料的缺乏，对其河型转化的控制因素的研究仍未得到系统的结论；本节将首先利用改进模型模拟游荡型小河的形成；然后通过调整其来水量、来沙量及边岸植被控制因素，进一步模拟该小河平面形态变化．通过模拟结果分析上述控制因素对游荡型小河转化的影响，并将结论与现有游荡河型转化理论进行对比．

图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 游荡型小河转化控制因素影响的数值模拟试验

基于上述产生的游荡型小河的地形，本小节通过概化数值试验，改变控制因素，模拟其河道进一步平面变形，试验分为 4组，试验条件及模拟结果见表2．计算过程中，水流时间步长为 6,s，泥沙时间步长为12,s，试验小河模拟历时600,d后地形高程见图4．

表2 数值试验条件及模拟结果Tab.2 Experimental conditions and simulation results

图4 试验小河600,d后平面形态Fig.4 Layout of experimental channel after 600,days

3.3 各控制因素对游荡河型转化影响分析

3.3.1 单一控制因素对游荡型小河平面形态变化的影响分析

图 4(a)～(c)分别为试验 A1～A3减小流量，减小来沙量及边岸植被覆盖的单因素改变下游荡型小河 600,d后的平面形态．试验 A1中减小流量相当于小水带大沙，此时的水流挟沙力小于其含沙量，因此河道上游的部分汊道将被淤死，出现了主槽；而在下游大量泥沙淤积于床面，使河道仍然散乱(图 4(a))；试验 A2减小来沙量，相当于大水带小沙，此时的河道挟沙力大于其含沙量，河道将冲深主槽，由于大水趋直，因此河道的下游出现了一段顺直河段，且汊道减少(图 4(b))；试验 A3中加入边岸植被影响，其边岸的抗冲性增强，限制河道进一步的展宽，由于边滩植草起到稳定边滩的作用，因此河型形态变化不大，一些汊道由于冲深而使水流归槽，从而其散乱程度有所降低(图4(c))．

选取试验 A1～A3组小河 6,000,m 处断面形态与游荡小河初始断面形态进行比较(图 5)：试验 A1的断面形态仍然为多汊道，但是其汊道数较初始为少；试验 A2的断面形态已经由原来的多汊道的“W”型转变为有主槽的“U”形；试验 A3小河由于边岸植被固岸的影响，因此展宽小于前两组，其形态与初始形态基本一致，主槽的刷深较初始地形加深，说明边岸植被有维持原有河道形态的作用．

图5 试验小河600 d后6,000 m断面形态对比Fig.5 Comparison of bed deformation at the 6,000,m sections of channel after 600 days

3.3.2 控制因素综合作用对游荡型小河平面形态变化的影响分析

图 4(d)为试验 A4控制变量(来水量，来沙量，边岸植被覆盖)综合作用下模拟小河 600,d后的平面形态．从表 2试验条件可知，随着来水来沙的减小，游荡型小河床面上的泥沙淤积减小，散乱程度减小，由于小水下水流归槽的作用，使得试验小河逐渐形成主槽；由于此时边岸植被固岸的共同作用下，小河的主槽能够得到进一步的维持，使其转化为弯曲小河的速度加快；当不考虑边岸植被影响只考虑来水来沙量减少时的游荡型小河转化为弯曲型小河需要的时间为 500,d，说明边岸植被起到稳定主槽作用，促进其游荡型小河向弯曲型小河转化．图6为试验A4模拟小河不同时刻的流场，可以看到来水来沙的减少，使水流逐渐归于主槽形成单一的河槽，从而促使其向弯曲型小河转化．

3.3.3 试验结论与现有游荡型河流转化理论对比

根据表 3[15]可知，当流量、来沙量减小，植被覆盖率增加的情况下，河型发展趋势为游荡向弯曲小河转化，通过数值试验 A1～A3可以得到：当仅改变单一控制因素时，游荡型小河的转化速度明显小于控制因素综合作用下的转化过程(试验 A4)；而通过试验A3对植被覆盖率对河型转化的影响模拟中，当其原河道较为散乱时，增加边岸植被可以稳定现有的河道，但是对其转化为弯曲河道的影响并不十分明显，这可能与模拟历时不够有关，但是在模拟结果的平面形态图4(c)中仍可以看到其汊道逐渐减少．

图6 试验A4小河不同时刻流场分布Fig.6 Temporal changes in flow field including velocity and surface elevation of test A4

表3 控制条件发生改变河型的发展方向[15]Tab.3 Transformation of channel with change of control factors[15]

4 结 论

(1) 考虑弯道二次流、河床结构对泥沙输运影响以及边岸植被影响的作用，改进原有模型，并提出简单易行的非黏性土崩岸模拟方法，成功模拟边滩形成及演变．

(2) 基于改进模型，建立概化游荡型小河转化控制因素影响的数值试验，首先模拟游荡型小河的形成；在此基础上，通过调整来水来沙量及边岸植被覆盖率等控制因素，分别模拟游荡小河进一步发展；通过模拟结果对控制游荡河型转化的因素进行分析得到：大水小沙、小水小沙、以及边岸植被覆盖率的增加都有使游荡型小河变为弯曲型小河的趋势；小水小沙以及加强边岸的植被覆盖率可以加速弯曲小河的形成；结论与现有河型转化理论吻合，说明利用数值模拟技术研究游荡河型转化及成因的可行性．

(3) 由于本模型考虑弯道二次流、河床结构及植被影响的公式与参数都是基于理论与试验相结合得到经验公式，因此泥沙的输运过程以及河道的横向崩岸过程将受到影响，从而导致河道平面形态模拟的不准确(如河道崩岸系数过大而导致河岸展宽过快造成流速降低，从而影响泥沙的输运及其进一步的形态变化)；同时，由于采用数组记忆边岸形态的变化过程，不能及时反映河岸崩塌过程对水流之间的相互关系．因此如何找到适应范围更广、理论上更合理的公式并改进数值计算方法，更为准确反映天然河道的演变过程将是未来河道演变模拟研究的重点．

［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］ 钟德钰，张洪武. 考虑环流横向输沙与河岸变形的平面二维扩展数学模型[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］ 周 刚，王 虹，邵学军，等. 河型转化机理及其数值模拟(Ⅰ)：模型建立[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］ 假冬冬．非均质河岸河道摆动的三维数值模拟[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］ 夏军强，王光谦.考虑河岸冲刷的弯曲河道水流及河床变形的数值模拟[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，Flokstra 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.