近半个世纪以来，随着计算流体力学理论及数值方法的飞速发展，冲积河流水沙数值模拟已成为研究河道中水流运动及泥沙输运不可或缺的手段，水沙数学模型也成为研究工程泥沙问题的重要工具之一^{[1, 2]}.根据边界上的水流、泥沙条件可以将水沙数学模型分为非恒定流模型和恒定流模型.天然河道中的水沙运动状态往往是非恒定的，在工程实践中为了节省计算工作量，通常将非恒定的来水来沙过程概化为梯级形式的来水来沙过程进行计算^{[3, 4]}，对于每一个梯级来说，流量为常数，水流近似为恒定流，与之对应的模型称为分级恒定流模型.

非恒定流模型与分级恒定流模型有其各自的优缺点，非恒定流模型的理论比较完善，因为其计算时间步长较短得到的结果更加准确，但相对耗时较长.而分级恒定流模型引入了“拟恒定流假设”，其适用范围有一定的局限性，适用于恒定或弱非恒定流的情况^[5]，而不能应用到像溃坝及冰湖溃决之类的强非恒定流情况^[6-10].但是相对非恒定流模型而言其计算时间步长较长，计算效率更高.

一维非恒定流模型有法国国家水利所的SEDICOUP模型、长江水利委员会长江科学院的HELIU-2模型、Change的FLUVIAL 11模型及Krishnappan的MOBED模型，一维分级恒定流模型有美国陆军工程兵团水文工程中心的HEC-6模型及Karim和Kennedy的IALLUVIAL模型等.Krishnappan^[11]用非恒定流模型MOBED及分级恒定流模型HEC-6预测加拿大Gardiner坝下游Saskatchewan河水流和泥沙输运及其河床高程的变化，并对两种模型的模拟结果进行了比较，发现非恒定流模型MOBED的模拟结果更接近实际情况.李义天^[12]等指出人为地将非恒定流概化为梯级式的恒定流求解在计算河段较短、河道的槽蓄量较小的情况下是基本正确的，但在计算河段较长、河道的槽蓄量较大的情况下往往不能正确模拟实际水沙运动过程，甚至无法直接应用.这只是一种直观的定性认识，其结论的正确性尚待验证.分级恒定流模型与非恒定流模型之间的定量差异，仍然缺乏足够的认识.本文通过非恒定流及分级恒定流数学模型来模拟在计算历时不断增加情况下概化长河段中的水沙运动过程，比较两种数学模型结果之间差异的变化.

1 数学模型 1.1 非恒定流模型

考虑一维矩形河道的挟沙水流过程，河床组成为无黏性均匀沙，泥沙的粒径为d.对水流、泥沙和河床变形建立质量和动量守恒方程，则完整的一维非恒定流水沙模型为

式中：U为守恒向量；F为通量向量；S_b为底坡源项；S_f为包括阻力项、含沙量梯度项以及水流与河床间泥沙交换项的源项；t为时间；x为平行于床面的空间坐标；h为水深；u为断面平均流速；c为体积含沙量；z_b为河床高程；g为重力加速度；S₀为阻力坡度；p为河床孔隙率；ρ_s为泥沙密度；ρ_W为水密度；ρ=ρ_W(1-c)+ρ_sc, 为浑水的密度；ρ₀=ρ_wp+ρ_s(1-p), 为河床饱和湿密度；E、D分别为泥沙的上扬和沉降通量.

对于上述非恒定流模型，本文采用有限体积法进行求解，首先通过算子分裂法将控制方程(1)离散为

式中：Δt为时间步长；Δx为空间步长；i为空间节点号；k为时间节点号；p为预报时刻节点号；F_i+1/2为x=x_i+1/2界面处的数值通量，本文采用具有TVD特性的SLIC算法进行求解^[13].

1.2 分级恒定流模型

在非恒定流模型的基础上将进口的来水来沙过程梯级化，假定一个梯级时间内的流量及输沙率恒定不变，并忽略式(1)中的∂U/∂t项，同时忽略河床冲淤变形对水流质量和动量改变的影响以及浑水动量方程右端的含沙量梯度项，就得到常见的分级恒定流模型.采用有限差分法离散可得

式中：q为单宽流量；z为水位；ψ₁、ψ₂分别为表征水流方程和泥沙方程的空间差分因子，通常要求0 < ψ_{1, 2} < 1.

1.3 封闭模式

应用曼宁糙率公式计算河床的阻力坡度：

式中：n为曼宁糙率系数；R=(B+2h)/(Bh), 为水力半径；B为河道宽度.

泥沙的上扬通量E及沉降通量D分别按下式计算：

式中：α为近底含沙量与平均含沙量之间的差异系数；ω为单颗粒泥沙在静止清水中的沉降速度，采用张瑞瑾沉速公式^[14]计算.

在本文算例中水流挟沙力的计算采用张瑞瑾公式，其显式表达式如下^[15]：

2 算例分析 2.1 算例介绍

计算河段长度选取长江中游河段从宜昌到九江的距离，取为900 km.将整个河段概化为矩形断面，河宽取为B，底坡为S_b，河床由无黏性均匀沙组成，泥沙粒径为d.在上游进口宜昌站给定来水来沙历时不断增加的情况下，来比较非恒定流模型与分级恒定流模型计算结果之间的差异.研究河道水流的实际情况往往涉及到很多不确定因素，本文研究的重点在于比较来水来沙过程梯级化后的分级恒定流模型与非恒定流模型之间的差异.作为初步研究，为了避免其他因素的干扰，暂时采用概化的地形条件.

在非恒定流模型计算中，假定洪水波来临之前河道水流近似为恒定流，采用一维恒定流水面线方程推求初始时刻各断面水面线的方法给定初始条件.上游边界条件为九江站的实测来水来沙过程，下游边界条件为九江站的水位流量关系，该水位流量关系曲线是根据九江站2004-2013年的水位流量数据拟合出来的，如图 1所示.进口断面的水沙过程选取宜昌站2013年的来水来沙过程作为代表年，为了研究在计算历时不断增加情况下两种模型计算结果之间差异的变化，将代表年2013年的水沙过程重复50次作为历时50 a的长系列年的水沙过程.

在分级恒定流模型计算中，上游边界条件为概化后梯级形式的来水来沙过程，下游边界条件也按水位流量关系给出.吴伟明等人^[16]在利用葛洲坝坝区1988-1991四年的水文泥沙资料对其所建立的模型进行验证时，将每年的来水来沙过程划分为80个梯级.本文参考其划分方法，在遵守来水来沙总量不变的前提下，用5 d一个梯级概化进口的来水来沙过程，来比较分级恒定流与非恒定流模型计算结果的差异.文中涉及到的分级恒定流模型若没特别说明，均是指按5 d一个梯级概化的分级恒定流模型.以历时1年为例，其梯级来水来沙过程如图 2所示.

在数值算例中，取河宽B=1 000 m，泥沙粒径d=0.062 5 mm，河道底坡S_b=5.5×10^-5，来水来沙总历时T=50 a.其他参数取值如下：泥沙的密度ρ_s=2 650 kg/m³，水的密度ρ_w=1 000 kg/m³，重力加速度g=9.8 m²/s，河床孔隙率p=0.4，糙率n=0.022，临界希尔兹数θ_c=0.04.河道底坡选取文献[17]中关于长江中下游宜昌大通河段的比降值，糙率值是根据初始时刻宜昌站的水位值来率定的.

2.2 冲淤量结果比较

表 1给出了通过非恒定流与分级恒定流模型计算河段冲淤量在不同时段的比较值，定性上来说，宜昌-城陵矶河段以及宜昌-九江总河段整体呈冲刷趋势，城陵矶-武汉以及武汉-九江河段先淤积后冲刷.当计算历时不断增加时，每隔10 a间河段的冲淤程度慢慢减弱，相应两种模型计算河段冲淤量之间的差异慢慢减小，但在时间上累计冲淤量之间的差异是不断变大的.从空间上看，当计算历时为50 a时，两种模型计算河段冲淤量在宜昌-城陵矶河段、宜昌-武汉河段及宜昌-九江总河段的绝对差异分别为0.45、0.72及0.80亿m³，相对差异分别为1.0%、1.4%及2.1%.随着河段距离的增加，两种模型计算河段冲淤量相差越来越大.

为了更加清楚地表现出两种模型计算河段冲淤量之间差异的变化，用分级恒定流模型计算河段冲淤量减去非恒定流模型计算河段冲淤量，得到两种模型计算河段冲淤量之间的差值，如图 3所示.图 3为分级恒定流模型中不同梯级概化方式计算河段冲淤量的差值比较图，从图中可以看出，在计算历时较短时，两种模型计算的总河段冲淤量之间的差异较小；随着计算历时的增加，差异变得明显.其中冲淤量差值为负值表示分级恒定流模型计算河段的冲刷量比非恒定流模型计算河段的冲刷量多.

以宜昌-九江总河段50 a的冲淤量为例，1 d一个梯级所对应的分级恒定流模型计算河段冲刷量要大于非恒定流模型的计算结果，冲淤量差值为-1.15亿m³.随着梯级时间间隔的变大，两种模型计算河段冲淤量的差异逐渐变小.当梯级时间间隔为15 d时，两种模型计算河段冲淤量的差值只有-0.04亿m³.之后随着梯级时间间隔的增加，河段冲淤量的差异又开始逐渐变大.当梯级时间间隔为30 d时，两种模型计算河段冲淤量差值为1.04亿m³，此时分级恒定流模型计算河段冲刷量要小于非恒定流模型的计算结果.

由于观测资料的限制，非恒定流模型给定的进口来水来沙过程与1 d一个梯级的分级恒定流模型的来水来沙过程相同.在非恒定流模型的计算过程中，考虑了水流、泥沙与河床的相互作用.当水流冲刷河床时，河床高程减小，过水面积增大，水流流速减小，从而减缓对河床的冲刷.但分级恒定流模型在一个时间步长内并没有考虑河床变形调整对水流的反作用，再加上流量差异对含沙量及河床变形的计算结果也有较大影响，导致1 d一个梯级的分级恒定流模型计算河段冲刷量要大于非恒定流模型的计算结果.对于分级恒定流模型来说，当进口来水来沙过程梯级时间间隔越长(即所分的梯级数越少)时，实际水沙变化过程越难体现出来，相应的梯级化过程越趋平缓，水流输运泥沙的能力越弱.所以梯级时间间隔长的分级恒定流模型要比梯级时间间隔短的分级恒定流模型计算的河床冲刷量小，1 d一个梯级的分级恒定流模型所计算的河段冲刷量是最大的.随着梯级时间间隔的增加，分级恒定流模型计算的河段冲刷量在逐渐减小的同时也在不断接近非恒定流模型的计算结果，此时两种模型计算河段冲淤量的差异是逐渐减小的；一旦分级恒定流模型计算的河段冲刷量小于非恒定流模型的计算结果，两种模型计算河段冲淤量的差异又开始逐渐增大.

实际工程中在对进口来水来沙过程进行梯级概化时，每一梯级所对应的时间间隔要适中，在本文中取15 d一个梯级较为合适，这样既能达到简化计算的目的，又可使分级恒定流模型计算河段冲淤量更加接近非恒定流模型的计算结果.

2.3 水位流量结果比较

图 4为非恒定流与分级恒定流模型在来水来沙历时1 a内计算结果在所选断面上水位及单宽流量随时间的变化过程图，从图中可以看出，两种模型计算结果在水位及单宽流量上均表现出一定的差异.非恒定流模型计算的断面水位变化过程与分级恒定流模型计算结果相比，存在相位落后的现象，且距进口断面越远这种差异就越明显.在x=650 km处水位差异最大为1.75 m、单宽流量差异最大为1.19 m²/s；在x=900 km处水位最大相差1.93 m、单宽流量最大相差3.54 m²/s.当计算河段较长时，洪水传播需要一定的时间，流量沿程会发生变化.而分级恒定流模型在某一时间梯度内，假定流量是沿程不变的，导致断面流量过程与非恒定流模型计算结果之间存在差异，因此分级恒定流模型无法预测某个时间梯度内洪水的传播过程.

2.4 河床纵剖面结果比较

图 5为非恒定流与分级恒定流模型计算河床纵剖面图，随着时间的增加，河段上游不断冲刷，下游不断淤积.为了更加直观地表现出两种模型计算河床高程之间的差异，用分级恒定流模型的计算结果减去非恒定流模型的计算结果，得到两种模型计算河床高程之间的差值ΔZ_b，如图 6所示.从图中可以发现，刚开始在下游出口处附近，分级恒定流模型计算的河床高程明显高于非恒定流模型计算河床高程值，在历时9.4 a时位于x=894.5 km处差异达到最大，最大相差3.35 m.这是因为两种模型计算的从上游携带的泥沙在下游堆积形成淤积体的位置之间存在一定的差异，这一点可以从图 5中的局部放大图中可以看出.随着时间的增加，淤积体位置不断下移直至泥沙从出口断面输运出去.之后两种模型计算的河床高程之间的差异开始在进口断面处体现出来，随后慢慢向下游传播，在历时34.3 a时进口断面处差异最大相差-0.32 m.

在上游进口附近，分级恒定流模型计算的河床高程值要低于非恒定流模型计算的河床高程值，下游出口附近分级恒定流模型计算的河床高程值要高于非恒定流模型计算的河床高程，这说明分级恒定流模型计算结果从上游往下游输送的泥沙更多.整体而言，河道上游发生冲刷，下游发生淤积，使得河床纵剖面向着减缓的趋势发展.而河床纵剖面减缓势必会降低河道水流输沙的能力，使得河道向着新的平衡状态转化.

2.5 计算时间比较

该数值算例中，非恒定流模型计算所需时间为84 695 s，分级恒定流模型计算所需时间为6 128 s，分级恒定流模型计算所需时间更短，其计算效率更高.在数学模型发展早期，由于计算资源的限制及计算机水平的影响，计算效率的问题成为科研工作者们优先考虑的因素，分级恒定流模型具有简化计算的优点就显得比较突出.但是随着计算机技术的发展以及程序结构的优化，非恒定流模型的计算效率也大幅度提高，科研工作者们在选取模型时更倾向于理论更加完善，更能反映实际水沙运动过程的数学模型，这个时候非恒定流模型就显示出其优越性.

3 结论

本文通过非恒定流及分级恒定流模型，模拟在计算历时不断增加情况下概化长河段中的水沙运动过程，并比较两种数学模型计算结果之间差异，所得结论如下：

1) 分级恒定流模型与非恒定流模型模拟的水位、流量等结果的变化过程在时间上存在相位差，分级恒定流模型无法描述水沙运动过程的时间效应，不能预测某一时间梯度内洪水的传播过程.

2) 分级恒定流模型计算总河段的冲淤量与非恒定流模型计算结果相比，在计算历时较短时差异较小.随着计算历时的增加，差异变得更加明显，当进口来水来沙历时长达50 a时，两种模型计算总河段冲淤量结果之间的相对差异为2.1%，且随着河段距离的增加，两种模型计算的河段冲淤量相差越来越大.在工程精度允许范围内，分级恒定流模型与非恒定流模型之间的差异可以忽略.分级恒定流模型可适用于短历时、短河段的情形；对于长历时、长河段水沙运动过程的模拟，推荐使用非恒定流模型.考虑到计算效率，如果误差在工程精度允许的范围内，则恒定流模型也是适用的.

3) 在概化梯级来水来沙过程时，当选取的梯级时间间隔过大或过小时，分级恒定流模型计算河段冲淤量与非恒定流模型计算结果相比均存在较大的差异，因此梯级时间间隔的选取要适中，本文中推荐选取15 d为一个梯级.

4) 虽然分级恒定流模型具有简化计算、提高计算效率的优点，但是随着计算机技术的发展及程序结构的优化，非恒定流模型的计算效率也不断提高.并且非恒定流模型的理论更加完善，更能反映实际的水沙运动过程，因此其应用前景更加广阔.