摘要：本文将正交曲线坐标系下的平面二维水沙数学模型和粘性河岸的冲刷模型结合，用于模拟弯道内的水流运动、悬移质泥沙的输移、河床的纵向及横向变形.用水槽试验资料验证了本文提出的水流模型，结果表明流速分布、水位等计算结果与实测值相当符合.应用建立的水沙数学模型以及河岸冲刷模型，模拟了一概化弯道在持续清水冲刷下的主流线位置、断面形态、主槽比降的变化过程，模拟结果符合弯道演变规律.

关键词：弯曲河道平面二维水沙数学模型河岸冲刷模型河床纵向及横向变形

弯曲河道是天然河道中常见的一种河流形态，多年来人们从各个方面对弯曲河道特有的水流运动规律、河床演变规律，进行了广泛的研究[1，2].如张红武[2]从模型试验和理论分析出发，较为深入地研究了弯道内的水面形态、环流运动以及纵向流速的沿程分布规律等.随着计算机性能的提高，以及数值计算方法的发展，很多研究者开始用数学方法模拟弯道内的水沙运动及河床变形.最早以engelund[3]、odgaard[4]等人提出的恒定水流模型为代表.近年来，部分研究者提出的三维数学模型，既能模拟弯道内的水沙运动，又能考虑床面的冲淤变化，较以前有了更大的进步[5，6].尽管三维数学模型可以模拟弯道内较为复杂的水流运动和河床变形问题，不过在实际工程中用的最广的还是平面二维模型[7，8].

通常现有这些模型只考虑河床的纵向变形，而很少注意弯道内的横向变形问题，也即河岸冲刷问题.本文针对当前在弯道水流和河床变形模型的研究状况，以前人的研究结果为基础，建立了正交曲线坐标下的平面二维水沙数学模型，用于模拟弯道内的水沙运动及河床纵向变形；在深入分析河岸冲刷机理的基础上，采用力学模型模拟粘性河岸冲刷、崩塌的过程.最后利用弯道模型的试验资料，验证了本模型计算的流场、岸边水位等与实测值符合较好.同时模拟了90°弯道在持续清水冲刷下的河床变形过程，对主流线、断面形态、主槽比降等计算结果随时间变化的分析表明，模拟结果符合弯道变形规律.

1数学模型

本文建立的模型具有如下特点：(1)为适应不规则的岸边界，建立正交曲线坐标下的平面二维水沙模型，并用madi方法[9]求解水流方程组.(2)改进以往对泥沙恢复饱和系数取为经验常数的做法，模型中采用了张红武[10]提出的不平衡输沙理论，用于研究悬移质泥沙的输移以及河床变形过程.(3)引入osman[11]提出的河岸冲刷模型，用于模拟粘性河岸的横向冲刷过程以及在重力作用下的崩塌过程.

1.1平面二维水沙数学模型

1.1.1水流基本方程

式中：εξ、εη分别表示ξ、η方向的泥沙扩散系数；s0k、sk、s*k、ωk分别为第k粒径组泥沙的侧向输入项、分组含沙量、分组挟沙力及有效沉速；α*、f1、k1分别为平衡含沙量分布系数，非饱和系数以及附加系数.

上述3个参数以及水流挟沙力的计算方法，详见文献[10].而分组挟沙力级配以及床沙级配的计算方法，见文献[12].

1.1.3河床纵向变形方程由悬移质泥沙不平衡输移引起的河床纵向变形方程为：

式中：ρ′为床沙干密度，n为非均匀泥沙的分组数.

1.2几个关键问题的处理

1.2.1正交曲线网格的生成设(x,y)为物理平面上的笛卡尔坐标系，(ξ,η)为相应计算平面上的变换坐标系，用以下方程可实现两个坐标系之间的变换，即：

求解上述方程组，即可得到物理平面(x,y)与计算平面(ξ,η)上对应点的关系.当网格正交时，β=0.控制函数的具体表达式，参见文献[13].

1.2.2初、边界条件对初始条件，一般给出水位、流速、含沙量等变量的初值.对进口边界、给定流量、沙量过程线，以及含沙量级配.对出口边界、给定水位过程线，并认为s/ξ=0.对岸边界，对水流可采用滑移或无滑移边界条件；对泥沙，可取s/=0(为岸边界法线方向).

在本文模拟的弯道中，横断面由主槽和滩地组成.假设滩地不过水，则岸边界是滩岸，同时认为滩岸可以冲刷.实际计算区域包括滩地和主槽，故采用“冻结法”技术[14]处理滩地上那些不过水的节点.一般岸坡坡顶和坡脚间的水平距离较小，可能远小于网格尺寸，在程序中不能分辨.为了能模拟河岸的冲刷过程，同时又不增加计算量，故本文在计算中用一数组记录岸边相邻两节点之间的实际地形.

1.3河岸冲刷模型现有的很多泥沙数学模型，很少考虑到河岸冲刷问题.即使那些考虑了河岸冲刷的模型，往往对弯道中河岸冲刷机理、冲刷速率，做了大量的假设和简化.如ikeda[15]认为弯道凹岸附近的纵向垂线平均流速大于弯道中心处的垂线平均流速时，弯道凹岸冲刷，反之则淤积.hasegawa[16]认为河岸冲刷速率与近岸的剩余流速成正比，根据泥沙连续方程，得出了适用于弯曲河道通用的河岸冲刷系数.不过，这类模型中存在一个共同缺点，即很少考虑到河岸组成物质、土体特性、几何形态等对河岸冲刷的影响.众所周知，粘性河岸和非粘性河岸的冲刷机理和冲刷速率是完全不同的.

为准确模拟河岸的横向冲刷过程，必须深入分析河岸冲刷机理.一般来讲，河岸冲刷后退是河岸土体和近岸水流相互作用的结果.除了岸边植被生长情况、河道内水位升降、渗流、管涌等因素影响河岸后退外，但主要是以下2种情况，是导致河岸后退的主要原因：(1)通过水流直接横向冲刷河岸导致河岸后退.近岸水流直接作用于河岸，冲动河岸边坡上水面以下的表层土体，并被水流带走，从而导致河岸后退.(2)通过河岸崩塌导致河岸后退.崩塌是河岸上的一部分土体在重力作用下，沿某一滑动面发生移动的过程.一般床面发生冲刷，导致河岸高度增加，或者水流淘刷河岸坡脚，使河岸坡度变陡，都会降低河岸的稳定性，当河岸的稳定性降低到一定程度后，河岸便会发生崩塌，导致岸边界后退.

图1粘性河岸冲刷的计算模式

本文主要考虑粘性河岸的冲刷后退过程，采用osman[11]提出的河岸冲刷模型(图1).该模型首先计算河岸横向冲刷距离，然后分析河岸是否会失稳、崩塌.在δt(sec)时间内，粘性河岸被水流横向冲刷后退的距离为：

式中：γs河岸土体的容重(kn/m3)；δb为δt时间内河岸因水流横向冲刷而后退的距离(m)；τ为作用在河岸上的水流切应力(n/m2)；τc为河岸土体的起动切应力(n/m2).cl为横向冲刷系数，取决于河岸土体的物理化学特性.

当由式(7)得河槽冲宽δb，用平面二维水沙模型算出河床冲深δz后，河岸高度增加，坡度变陡，稳定性降低.根据土力学中的边坡稳定性关系，采用若干假定，可得到河岸发生初次崩塌时的临界条件.若河岸已发生初次崩塌，则假定以后的河岸崩塌方式为平行后退，即崩塌后的边坡角度恒为β，仍可用土力学的方法判断是否会发生二次崩塌.

1.4计算方法和求解过程首先，采用madi法[9]计算流场.这种方法采用非交错网格，将水位、流速等变量均布置在同一网格点上，对水流连续方程和动量方程均不按方向剖开，由此将基本运动方程离散而形成新的差分代数方程组，并建立一种新的解法.具体求解过程如下：对式(1)进行差分离散，时间导数项采用前差表示，空间导数项采用中心差分表示，将时间步长δt分成两部分.在前半时间步长内，利用连续方程(1)、ξ方向动量方程(2)，沿ξ方向采用隐式离散，对η方向采用显式离散，求解得zn+1/2,un+1/2，再利用η方向动量方程(3)，显式求解得vn+1/2.在后半时间步长内，利用式(1)、(3)，沿η方向采用隐式离散，对ξ方向采用显式离散，求解得zn+1，vn+1，再利用式(2)，显式求解得un+1.由于要进行长时间的河床变形计算，计算量特别大，在实际计算中通常用伪恒定方法计算出恒定流场.然后，采用显隐混合格式计算悬移质含沙量分布.先按不同方向ξ、η对式(4)进行算子分裂，在前半步长内，对ξ方向的分步方程用指数格式显式离散求得sn+1/2，在后半步长内，对η方向的分步方程用crank-nicholson型格式隐式离散求得sn+1.以时间为迭代参数，计算出恒定的浓度场.

最后，显式求解式(5)，可得出时段末的床面高程.在已知河床变形的基础上，采用河岸冲刷模型，模拟弯曲河道河岸冲刷后退的过程，具体计算方法见文献［12］.此外，本文在计算中还采用以下假定：直接从河岸冲刷下来的和河岸崩塌后产生的泥沙，全部转化为悬移质泥沙，并作为下一个时段的侧向输沙量.

2模型的验证

为了准确模拟弯道的河床变形过程，本文首先对水流模型进行验证.验证的资料来自罗索夫斯基[1]的180°弯道的水槽试验资料.该水槽由一6m长的顺直进口段、180°的弯道段以及3m长的出口段组成.水槽宽度为0.80m，过水断面为矩形.弯道外半径为1.2m，内半径为0.4m.计算中的水流条件如下：流量为12.3×10-3m3/s，进口断面选在水槽进口以下2m处；出口断面选在距水槽出口0.1m处，断面水深为0.054m.水槽糙率取0.012，计算区域内共划分网格数为95×11.

图2弯道凹岸与凸岸岸边水位的沿程分布

图3纵向垂线平均流速分布

图2给出了弯道中凹岸和凸岸水位沿程变化情况.从图中可以看出，在弯道进口段，凹岸和凸岸的水位基本相同；在弯道段，凹岸水位明显高于凸岸；在弯道出口段，两岸水位基本相同.尽管计算值略小于实测值，但两者的变化趋势相当符合，当水流由顺直河段进入弯道后，由于受到离心力的作用，使弯道凹岸一侧的水位恒高于凸岸一侧，最大水位差一般出现在弯道顶点附近，而向上下游两个方向逐渐减少.由于未能考虑弯道内横向环流对水流动量方程的影响，从而导致水位计算值总体偏小.

图3给出了纵向垂线平均流速的沿程分布情况.从图可知，实测值与计算值符合较好，尤其在弯道的进口段附近.从图中还可以看出，在弯道的进口段，最大纵向平均流速的位置(主流线)紧靠凸岸；在弯道段，主流线仍紧靠凸岸；在弯道出口段附近，主流线逐渐向凹岸转移.出弯后的水流，在相当长的距离内，最大流速仍靠近外侧(凹岸).