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

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

弯曲河道是天然河道中常见的一种河流形态，多年来人们从各个方面对弯曲河道特有的水流运动规律、河床演变规律，进行了广泛的研究[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水流基本方程

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

上述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]处理滩地上那些不过水的节点.一般岸坡坡顶和坡脚间的水平距离较小，可能远小于网格尺寸，在程序中不能分辨.为了能模拟河岸的冲刷过程，同时又不增加计算量，故本文在计算中用一数组记录岸边相邻两节点之间的实际地形.

3计算结果与分析

由于目前缺少由悬移质泥沙不平衡输移的实测资料，故本文采用一概化的90°弯道，研究其在持续清水冲刷、两岸边界可动的条件下，弯道内流速、断面形态，以及比降的变化过程.对模拟结果的分析表明：本文所建立的平面二维水沙数学模型以及河岸冲刷模型，能较好的反映弯曲河道的演变过程.

3.1河道概况假定有一概化弯道，包括顺直进口段(1000m)、90°弯道段(2100m)以及顺直出口段(2000m).河床初始纵比降为0.0002，初始糙率为0.02，床沙平均粒径为0.087mm.河道主槽宽度约450m，两岸滩地总宽约250m；初始河岸坡度和高度分别为60°、5.0m(见图5).河岸土体特性为：容重γs=18kn/m3、内摩擦角φ=14°、凝聚力c=20kn/m2、τc=1.2n/m2.模拟河段共划分网格为50×28，每个计算时段恒为2d.假设该弯道进口处的悬移质含沙量很小，计算中取0.10kg/m3，而流量始终为3000.0m3/s；弯道出口断面初始水位为104.0m，在水流的持续冲刷下，水位不断下降，在计算中假定水位下降幅度为河床冲刷厚度的0.2倍.

3.2流速分布图4为初始时刻的流速分布图.在初始时刻，主流线在进口处靠近凹岸，当水流进入弯道后，由于离心力的作用，主流线逐渐贴近弯道凸岸，在弯道出口处，主流线又逐渐偏离凸岸，而向凹岸靠近.当该弯道经过30d的持续清水冲刷后，平均水深3.9m增大4.7m，平均流速由1.6m/s减小1.2m/s.不仅弯道内的水流条件发生变化，而且断面形态、平面形态也发生变化，从而使弯道内的主流线位置也发生了变化.从图5中可以看出，在弯道进口段，主流线基本上贴近凹岸；在弯道段，主流线仍贴近凸岸；在弯道出口段，主流线逐渐向河道中心位置转移.弯道内水深增加，流速减小，从而导致主流顶冲凹岸的位置逐渐上移.主流线顶冲凹岸的部位，符合小水上提，大水下挫的弯道水流运动规律.此外，从图5中还可看出，由于河岸冲刷、崩塌后退，原先的那些不过水的滩地位置，逐渐开始过水.

由此可见，弯道内垂线平均流速的分布，不仅与进口处的流速分布、出口断面的水位有关，而且还与弯道的断面形态、平面形态等因素有关.

图4初始流速分布

图530d后的流速分布

图6断面形态变化过程

图7主槽平均高程变化过程

3.3断面形态变化图5给出了弯道段第21号断面的断面形态随时间变化过程.从变化图中可以看出，由于两岸滩地可以冲刷，主槽在不断冲刷下切的同时，河床逐渐展宽.由于在弯道段，主流线的位置一直靠近凸岸附近，导致凸岸附近的流速大于凹岸附近流速，因而凸岸的滩地冲刷、崩塌后退的距离(28.6m)略大于凹岸滩地后退距离(28.4m).而且随着河床逐渐展宽，滩地冲刷、后退的速率逐渐减小.另外，30d内河床纵向总冲深为1.37m，每隔10d的冲刷深度分别为0.61m、0.43m、0.33m.由此可见，在进口水沙条件不变的条件下，河床的纵向冲刷速率随时间增加而逐渐减小.

弯道进口段和出口段的断面形态变化，也有类似的变化规律.不过，应当指出，这里未出现一般弯道演变中的凹岸冲刷，凸岸淤积的现象，主要由以下两方面的原因：一是河床受清水冲刷，凸岸附近的含沙量(平均为0.27kg/m3)远小于当地的水流挟沙力(平均为3.2kg/m3)，因此发生冲刷.二是由于水流中挟带的悬移质含沙量本身就很少，由横向环流从凹岸带到凸岸的悬沙数量也少，加上凸岸没有高程较低的边滩，因而凸岸不发生淤积.

3.4主槽比降的变化图6给出了河岸可冲刷情况下，主槽平均高程随时间的变化过程.在清水的持续冲刷下，整个河床都有不同程度的冲深，30d内，进口处冲深较多(1.3m)，出口处冲深相对较小(1.0m).主河槽纵比降初始值为2×10-4，河床冲刷10d、20d、30d后的比降变为1.82×10-4、1.63×10-4、1.59×10-4.随着河床的逐渐冲深和横向展宽，比降变小的速度降低.本文还计算了河岸不可冲刷情况下，主槽平均高程随时间的变化过程.在河岸不可冲刷时，河床冲深明显加大.在30d内进口处冲深约1.5m，出口处约1.2m.河床纵比降由原来的2×10-4下降为1.4×10-4，下降幅度为30.

从上述分析来看，在河岸可冲刷的情况下，由两个重要原因可以减缓河床冲刷：一是河槽宽度增大，导致过水面积增大，流速较小，水流挟沙力减小，从而减少对床面的冲刷.二是从河岸冲刷和崩塌下来的泥沙数量部分满足了悬移质水流挟沙力的要求，因而可减少从床面冲起的泥沙数量.此外，还有河床冲刷过程中产生的床沙粗化现象，可引起河道阻力增加，也可减缓对河道的冲刷.

4结论

本文在分析前人对弯道内水流运动和河床变形的基础上，针对当前很多模型中较少考虑河岸冲刷问题的状况，提出了正交曲线坐标系下的平面二维水沙数学模型，结合osman提出的粘性河岸冲刷模型，用于模拟弯道纵向及横向的河床变形.通过对一水槽弯道模型和一概化弯道的模拟，本文得到以下结论：

(1)采用曲线坐标下的水流控制方程，能有效的解决不规则岸边界问题；用madi法求解方程，不仅能减少计算时的存储量，而且减少计算量.通过对水槽模型试验资料的对比，本文提出的水流模型的计算结果具有一定的精度和稳定性.

(2)采用非均匀悬移质泥沙输移方程，引入张红武提出的非平衡输沙模式，克服了以往对方程中的泥沙恢复饱和系数任意调整的局面.

(3)平面二维水沙数学模型与河岸冲刷模型相结合，能更好的模拟河道的河床变形过程.通过对一概化弯道的模拟，表明弯道内主流线变化不仅与进出口水沙条件有关，而且还与断面形态有关.而断面形态的变化与河岸抗冲性能有关.对于弯道在持续清水冲刷下，当河岸可冲刷时的主槽比降变化比河岸不可冲刷时要慢.

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