第二节 泥沙数学模型

泥沙数学模型计算是三峡工程泥沙问题研究的主要方法之一，在三峡工程论证和设计阶段为工程的可行性论证和运行调度方案的确定提供了重要依据。三峡工程泥沙数学模型计算研究大体分为两个阶段：1958—1971年阶段和1972年以后阶段。

1958—1971年阶段长江科学院主要采用由输沙平衡原理建立的有限差分法即平衡输沙法，以及平衡坡降法即三角洲淤积计算法等计算方法，进行三峡水库泥沙淤积量和淤积年限计算。

1972年以后，长江科学院韩其为、黄煜龄、欧阳渊等在武汉大学数学系杨曙光老师协助下以悬移质不平衡输沙研究成果为基础，开始研究水库泥沙冲淤计算方法和研制电子计算机计算软件，经过一年多的摸索，编制调试出一套水库泥沙冲淤计算通用程序，并利用丹江口水库1967—1969年水库淤积资料和荆江严家台放淤区1966年淤积资料进行了验证^［4］。以后软件研发人员对上述水库泥沙淤积计算通用程序作了进一步改进和计算机调试，成为“HELIU-1”软件，可用于水库泥沙淤积计算^［5］。应用该软件进行了三峡工程正常蓄水位150m、160m、170m、175m、180m方案的水库泥沙淤积计算。为配合三峡工程可行性重新论证阶段有关坝下游河道冲刷问题研究，1985—1990年长江科学院黄煜龄、黄悦等研制了河道冲淤计算程序，即“HELIU-2”软件，该程序可用于河道冲淤计算和水库泥沙淤积计算^［6-8］。

1986年三峡工程可行性重新论证阶段，中国水利水电科学研究院韩其为等在以往数学模型研究的基础上，进一步发展为M1-NENUS-2模型和M1-NENUS-3模型，用于三峡水库泥沙淤积和坝下游河道冲刷计算^［9-11］。对于前者，1986年三峡工程泥沙论证专家组林秉南院士认为^［12］：①该模型适用于小含沙量的流动，长江含沙量在20kg/m³以下，可视做小含沙量；②该模型适用于细颗粒泥沙，川江悬移质中值粒径为0.034～0.036mm，可看做细颗粒泥沙；③该模型有一定的经验性，如恢复饱和系数α、出现平衡糙率时的断面积a_k、冲淤在断面中的分配或分布形式以及床面形态糙率等取值或假定来自以往的经验和试验及实测资料，从已有验证资料来看，采用同一组系数已在不同情况下取得较好验证，现阶段应用这些系数和假定应是可行的，在可行性研究阶段，应用该数学模型对长河段的冲淤作出的估计也应是可行的。

随着三峡工程技术设计及工程建设的需要，在水库泥沙淤积与水库优化调度等课题研究中，促进一维数学模型由恒定水沙模型发展为非恒定水沙模型，如清华大学水利水电工程系周建军、中国水利水电科学研究院毛继新及长江科学院黄仁勇等，于本世纪初先后提出了适用于三峡水库泥沙冲淤计算的一维非恒定水流泥沙数学模型。

在三峡工程泥沙问题研究中，南京水利科学研究院窦国仁、武汉大学李义天、长江科学院范北林、中国水利水电科学研究院方春明等先后建立的二维泥沙数学模型共10多个，多用于水库变动回水区冲淤计算和坝下游河道冲刷计算。这些泥沙数学模型采用的基本原理相同，方程组基本一致，由于研究目的不同或研究对象的差别较大，各模型的研制者针对不同情况，采用的计算方法也不尽一样。

21世纪初泥沙数学模型研究发展较快，不仅二维泥沙数学模型在计算方法上呈现多样化，三维泥沙数学模型也开始应用于三峡工程泥沙问题研究，如南京水利科学研究院陆永军等建立了以紊流随机理论为基础的三维泥沙数学模型。

三峡工程可行性重新论证和设计阶段，一维、二维泥沙数学模型研究成果较多，见表2-4。限于篇幅，本节仅就有代表性的一维泥沙数学模型作简要介绍。

一、一维恒定非均匀不平衡输沙数学模型（M1-NENUS-3）

M1-NENUS-3模型是中国水利水电科学研究院韩其为等在M1-NENUS-2模型基础上建立并完善的一维恒定非均匀不平衡输沙微冲微淤数学模型。该模型能较好地模拟河道冲淤和水库淤积，曾在“七五”、“八五”国家重点科技攻关项目、“九五”三峡工程泥沙问题研究中得到较多的应用^［7-9］。

（一）基本方程

1.水流运动方程

式中：H为水位，m；n为糙率；Q为流量，m³/s；A为过水面积，m²；B为水面宽度，m；Δx为断面间距，m；g为重力加速度；i为时段编号；j为断面编号（自上向下）。

2.水流连续方程

（1）流量为常数时

（2）流量近似为沿程线性变化时

式中：下标“0”、“1”分别表示计算河段的进、出口断面；L为河段总长度，m；其余符号意义同前。

3.悬移质含沙量方程

其中

以上式中：S为含沙量，kg/m³；S^*为挟沙能力，kg/m³；P为悬移质级配；P^*为挟沙能力级配；ω^*为挟沙能力级配相应的沉速，m/s；α为恢复饱和系数，淤积时取0.25，冲刷时取1.0；脚标“L”表示泥沙的组数；ω（L）为第L组泥沙的沉速，m/s；其余符号意义同前。

4.河床变形方程（泥沙连续方程）

式中：Δa_i，j为j断面冲淤面积，m²；Δt_i为i时段时间步长，s；为淤积物干容重，kg/m³，初期淤积物干容重由淤积物级配决定，计算年限较长取密实后干容重。

（二）辅助方程

1.挟沙能力及挟沙能力级配方程

挟沙能力级配指与挟沙能力相应的级配，它等于在同样水流和床沙条件下，输沙达到平衡时的悬移质级配^［37-39］，其方程可按三种输沙状态来确定。

（1）明显淤积。所谓明显淤积是指各组粒径泥沙都发生一定程度淤积（至少不发生冲刷）。这种淤积在床面淤积速度较快时出现。从瞬时情况看，在淤积过程中床面泥沙中虽有被冲起但冲起后又被淤下，从累计效果看，床面泥沙不会被冲起，即此时的床面变形和悬移质运动与原河床无关。在明显淤积条件下，可从理论上证明^［37］

（2）明显冲刷。明显冲刷是指各组粒径泥沙都发生一定程度的冲刷（至少不发生淤积）。这种冲刷在床面冲刷速度较快时出现。从瞬时情况看，在冲刷过程中悬移质泥沙虽有被淤下的，但淤下后又被冲起，从累积效果看，悬移质不会被淤下，即含沙量和悬移质级配的沿程变化单纯由沿程冲起泥沙数量和级配而引起。明显冲刷条件下，有近似关系^{［37，38］}

明显冲刷和明显淤积状态下的挟沙能力级配等于或约等于悬移质级配，统称为明显冲淤。以这种挟沙能力级配所建立的模型相应称之为明显冲淤模型，如M1-NENUS-2模型。其挟沙能力公式为

而

上二式中：K₀为挟沙能力系数，需根据实际河段而定；m为指数，m=0.92。

（3）微冲微淤。微冲微淤是指各组泥沙的冲淤性质可能不一样，有几组泥沙被冲起，而另有几组可能发生淤积。微冲微淤时挟沙能力级配与悬移质级配不完全一致，挟沙能力级配不但跟悬移质级配有关，还跟床沙级配有关。根据文献［37、38］给出微冲微淤挟沙能力级配如下：

其中

以上式中：脚标带“0”的为相应参数在短河段进口断面的值；S^*（L）为第L组泥沙的均匀沙挟沙能力，即用均匀沙沉速ω（L）代入式（2-9）、（2-10）即可；P_4，L，0为进口断面悬移质级配；k为悬移质中的细颗粒组数；N为悬移质总组数。

以上式中：P_1，L，1为床沙级配；n为床沙中可悬浮泥沙组数；M为床沙总组数；为与该组床沙级配P_{4，L，1，1}相应的挟沙能力级配；S^*（）为其相应的挟沙能力。

取从冲刷开始（λ^*=0）至冲刷后（λ^*）与粗化床沙P_{4，L，1，1}相应的平均挟沙能力级配，即

此处

而ω_m由

试算确定。

从式（2-11）可知，微冲微淤挟沙能力级配由三部分组成：第一部分为悬移质中的细颗粒，这些来沙的累积效果是不淤的，因此它的挟沙能力就是这部分来沙；第二部分为悬移质中的较粗颗粒，这一部分泥沙转成床沙后，只有部分能转成挟沙能力；第三部分为从床沙中冲起来的部分。微冲微淤挟沙能力公式为^［38］

式中符号意义同前。

水库下游冲刷过程一般都是细沙冲、粗沙淤的分选过程。当水流强度不是很大时，各组泥沙有冲有淤，冲淤性质可能不完全一样。因此，从物理模式来看，微冲微淤模型更能反映实际情况。

2.悬移质级配变化

（1）明显淤积。明显淤积时悬移质级配为

其中

以上式中：θ为修正指数，湖泊型水库取θ=1/2，河道型水库取θ=3/4；λ为淤积百分数，由式（2-16）输沙率定义；ω_m，i，j为中值沉速，由式（2-27）确定。

（2）明显冲刷。明显冲刷时的悬移质级配为^［39］

其中

式中：P_{1，L，i-1，j}为该计算时段冲刷开始时的床沙级配；ω_m，i，j由确定；λ^*为冲刷百分数，表示泥沙冲刷分选程度；Δh′为虚冲刷厚度，m，表征单位面积上冲刷泥沙重量（t/m²）；Δh₀为参加冲刷分选的厚度，m，表征意义同上；B_k为稳定河宽，m；其余符号意义同前。

（3）微冲微淤。微冲微淤时悬移质级配的确定是先由式（2-4）求出分组含沙量，再求悬移质级配。

（三）微冲微淤模型的使用范围

微冲微淤是指一次冲刷过程的冲刷幅度较小，各组泥沙有冲有淤，这与明显冲淤的性质是完全不同的。微冲微淤的使用范围是式（2-11）和式（2-24）必须满足

当不满足式（2-31）时，模型自动转入M1-NENUS-2明显冲淤模型计算。

（四）模型验证

采用丹江口水库1968—1986年淤积资料对M1-NENUS-2模型进行了验证，验证计算河段从白河水文站至大坝，包括整个汉江干流库区，长约200km。采用1981—1987年长江中游宜昌至城陵矶河段冲淤资料对M1-NE NUS-3模型进行了验证。

二、一维全沙河床冲淤数学模型（HELIU-2）

长江科学院黄煜龄、黄悦等在“七五”期间提出的一维全沙河床冲淤数学模型，经“八五”期间逐步加以改进与完善，主要应用于河道长河段泥沙冲淤问题研究，20世纪90年代中后期该模型进一步拓展应用于大中型水库泥沙淤积问题研究。该模型能较好地模拟悬移质和推移质运动条件下长河段的河床冲淤和水库淤积规律，曾在“七五”、“八五”国家重点科技攻关项目，“九五”三峡工程泥沙问题课题研究中得到较好的应用^［10-12］。

（一）一维全沙数学模型基本方程

1.水流运动方程和连续方程

2.悬移质泥沙连续方程

3.水流挟沙力方程

4.推移质输沙率方程

5.悬移质河床变形方程

6.推移质河床变形方程

以上式中：Z为水位；Q为流量；J_f为能坡；U为流速；A为过水面积；g为重力加速度；S和S_*分别为断面平均含沙量及挟沙力；G_b为推移质输沙率；γ′_s为淤积物干容重；ΔA₁、ΔA₂分别为悬移质和推移质冲淤面积；B为水面宽；x为沿程距离；t为时间；ω为泥沙颗粒静水沉速；脚标“i”为第i粒径组；d为粒径；α为恢复饱和系数（选用0.25）。

（二）基本方程组简化

由于所研究的问题是长时段、长河段内发生的冲淤，在实际计算时对基本方程组进行了简化，将整个计算时段划分成若干个小的计算时段，将长河段划分为若干个短河段，且在计算时段内和河段内除ΔA以外，其他因子不变，即按恒定流考虑；而在不同时段不同河段各因子可以不同。经简化和推导，可得到应用方程组。

1.水面线计算式

2.悬移质含沙量变化方程

分组挟沙力系数K_i，采用窦国仁公式计算^［24］：

悬移质级配P_i，由下式计算：

3.悬移质运动引起的河床变形

4.推移质输沙率

推移质输沙率用长江科学院提出的输沙经验曲线（图2-6），输沙曲线的关系式为

5.起动流速

起动流速采用张瑞瑾公式计算^［41］：

6.推移质运动引起的河床变形

7.河床总变形

8.床沙级配计算

设G_pi为分组冲淤量（kg），冲为“-”，淤为“+”；G_ci为河床分组可冲量（kg），“+”；G_ti为河床分组剩余量（kg），“+”；若G_pi≥0，则

G_ti=G_pi+G_ci

若G_pi＜0，G_ci≥|G_pi|，则

G_ti=G_ci+G_pi

若G_ci＜|G_pi|，则

G_ti=0

故床沙级配为

9.横断面形态修改

按沿湿周等厚变形修改。淤积时，全断面淤积；冲刷时，冲槽不冲滩。

10.宽断面的处理

宽断面主要由滩和深槽两部分组成。中低水时，边滩部分基本不过流；高水时，主流走深槽，滩上的水深小，主要起蓄水作用，是泥沙淤积的部位。因此在河床变形计算时，对宽断面分滩槽两部分作简单处理，主槽河宽取相当造床流量时河宽，在验证计算时以各分河段的最小河宽与划分的计算河段内的平均河宽之间调整，作为滩槽分界点。

以上式中：Δt为时段；Δx为两断面间距；S_i、S_*i分别为分组含沙量及挟沙力；S_*m为断面总的挟沙力；q为单宽流量；L为河段长度；ω_m为非均匀沙平均沉速；ω_i为第i组泥沙颗粒沉速；K和m分别为挟沙力系数和指数；β为指数，取1/6；v_d为近床面流速；v_c为床沙起动流速；d₅₀为中值粒径；h为平均水深；G_s为悬移质输沙率；g_b为推移质单宽输沙率；G_b为推移质输沙率；脚标“0”代表已知断面；其余符号意义同前。

（三）数学模型的解法

在含沙量不大的河流，假设泥沙颗粒运动不受相邻颗粒运动的影响，可把混合沙的每一组粒径单独计算，然后把各粒径组的部分成果加起来，就得到总的结果。

基于上述假设，求解时先求解各粒径组的含沙量、冲淤量及所引起的河床变形，然后将各粒径组的计算结果合起来就得到河床冲淤变形的总结果。由于按粒径分组计算冲淤量，在同样的边界与水力因素的条件下，能自动进行悬移质中的粗颗粒与河床中的细颗粒交换，达到“淤粗悬细”的结果，在冲刷过程中自行完成粗化过程，形成粗化保护层。

求解方程组时采用非耦合解，每个计算时段分三步计算：第一步推求水面线，算出各断面的水力要素；第二步求各河段各组泥沙（包括推移质和悬移质）的冲淤量；第三步修改横断面形态。

长江科学院用上述方法研制了河道冲淤变形计算软件（HELIU-2程序）。程序采用模块化编制，流程见图2-7。

（四）数学模型的验证

采用长江宜昌至大通长1125km河段1980—1987年河道冲淤资料、汉江丹江口水库汉江库区1968—1985年淤积资料以及坝下游襄阳至仙桃长320km河段的冲淤资料对HELIU-2模型进行了验证。

三、一维非恒定水流泥沙数学模型（周建军）

清华大学周建军建立的一维非恒定水流泥沙数学模型是依据不恒定水流数学模型和不恒定、不均匀的悬移质泥沙数学模型展开的。在“十五”三峡工程泥沙问题研究中应用于三峡水库淤积计算^［16］。

（一）数学模型基本方程

水流连续方程

水流动量方程

分组不平衡输沙方程

河床冲淤变形方程

其中，Φ_k=P_bkφ_k+（1-P_bk）min（S_k，φ_k）是第k组泥沙的冲刷函数；φ_k=P_kφ是分组挟沙能力；

是混合沙的挟沙能力；；K_s为混合沙挟沙力系数；P_k是分组挟沙能力级配，根据李义天的方法计算：

其中，，，；P_bk是分组床沙级配。在式（2-50）中，G为推移质输沙率。床沙组成采用分层储存级配、表层参加交换的计算模式。表层床沙级配变化按下面公式计算：

其中，Δ_b=2.5δ_e；；；δ_k是各个分组在时段的实际冲淤厚度；是时段初的床沙级配。

（二）综合恢复饱和系数及断面冲淤修正方法

1.综合恢复饱和系数方程

式（2-49）中是综合恢复饱和系数，与河道断面形态、冲淤状态及水沙运动有关，据作者研究：

其中

式中：β由方程

的第一个非零正根确定；R=6ω/κu_*；κ是Karmann常数；u_*是摩阻系数；α^*是相当于均匀水流中的泥沙恢复饱和系数，它不但决定于水流泥沙参数R，而且与冲淤状态关系很大。其变化情况如图2-8所示。

式（2-54）确定的综合恢复饱和系数反映了水流泥沙条件及河道断面诸多因素的作用，是适合于一般情况的参数。在矩形渠道中，α^s=α^*，恢复饱和系数都是大于1.0的。而且，冲刷情况下的恢复饱和系数大于淤积时，公式中η=h/H，h是垂线水深，是断面平均水深；v是反映断面水流速度分布的参数，即横向流速分布与水深分布间的幂指数关系：

其中

式中：B是河道水面宽度。断面形态对综合恢复饱和系数的影响与断面形状和水位等有关系。

2.冲淤断面修正方法

严格来讲，一维数学模型不能准确计算泥沙冲淤在横断面上的分布。但是，河道断面的冲淤分配是模型必需的一个环节，它在很大程度上影响冲淤计算的成果。迄今所有的一维数学模型都采用经验方法修改河道断面，比如冲刷和淤积按湿周分布，或者淤积沿断面均匀分布，冲刷只影响主槽等。经验方法具有一定的依据，但也给数学模型带来了较大的任意性。特别是在河道淤积和冲刷交替的情况下，不能塑造合理的平衡断面，甚至将影响数学模型计算的结果。该模型依据流管积分和恢复饱和系数的理论体系，提出了冲淤断面修正的理论公式：

式中：δZ_b和δA_s分别为垂线冲淤厚度和断面冲淤面积。由此可得断面冲淤分配系数为

式（2-59）给出了基本符合河道冲淤经验的断面形态，即淤积沿河宽分布相对比较均匀，而冲刷主要集中于深槽（见图2-9）。

（三）计算方法

该模型不恒定水流的计算方法采用了Preissmann四点隐格式。该方法具有很好的稳定性和精度。与常用的推求水面线方法相比，该方法可以很好地适应缓流、急流及其过渡计算，对河道中峡谷段可能出现的急流能很好地模拟。

关于不平衡输沙方程，在不恒定流条件下，采用迎风差分格式进行计算。

（四）基本参数的确定方法

在一维数学模型中，基本参数包括河道糙率及其变化规律、泥沙沉降速度、挟沙能力系数和指数、泥沙淤积物的干容重等。

蓄水后的天然河道糙率可采用2003年三峡水库蓄水到135m后的沿库水位资料率定。水库淤积后，完全覆盖的床面糙率根据冲刷河道的糙率资料计算（采用韩其为在论证期间计算采用的荆江河道的糙率资料）。关于天然河道向冲积河道的糙率过渡，采用下列插值方法：假设天然河道床面突体高度δ_b、天然河道糙率n₀和冲积床面糙率n_f，根据断面上点的淤积厚度，按线性插值可以给出该点不同时候的糙率：

当δZ_b>δ_b时，n=n_f。而断面平均糙率为

式中：dp为断面湿周微分；χ=∫dp为断面湿周；计算中取δ_b=10m。

采用该方法计算断面综合糙率与δ_b的选取有关，但是三峡水库淤积量较大，其影响持续时间不长。这与黄煜龄和韩其为用淤积面积为参数直接插值相比，本方法给出的过渡期糙率相对要小些，这是符合实际情况的。实际上，当床面被泥沙覆盖但流速还没有大到可以在河床形成沙波等形态时，水库的糙率应该比冲积河道的糙率更小。

泥沙沉降速度计算，采用水利部规范指定的沉降速度公式。计算沉降速度时，水温按23℃计，泥沙比重γ_s=2.65。挟沙能力按武水院公式计算^［41］，其中挟沙能力系数和指数按韩其为的方法取为K_s=0.246，m=0.92，经过2003年资料检验和调整后，也采用了较小的挟沙力系数计算；泥沙水下干容重按1325kg/m³计算；式（2-57）中，断面水流速度分布参数v，根据Manning公式的结构形式取值为0.667。

（五）模型验证

选用三峡水库2003年蓄水运行一年的水库淤积资料对本模型进行验证计算。计算范围：干流上游从距大坝756km（朱杨溪以上10km）开始到三峡大坝。

四、一维非恒定水流泥沙数学模型（毛继新）

中国水利水电科学研究院毛继新提出的一维非恒定水流泥沙数学模型建立于三峡工程试验性蓄水阶段。该模型应用于“十一五”三峡工程泥沙问题研究中的水库淤积计算^［17］。

（一）数学模型基本方程

水流连续方程

水流动量方程

河床变形方程

泥沙连续方程

水流挟沙能力

以上式中：Q为流量；A为断面面积；B为断面宽度；Z为水位；K为流量模数；S为断面平均含沙量；S^*为断面平均挟沙能力；g为重力加速度；α为泥沙恢复饱和系数；ω为泥沙颗粒沉速；ΔA为断面冲淤面积；γ_s为淤积物干容重；k为挟沙能力系数。

（二）方程求解

一维非恒定流泥沙数学模型的计算采用非耦合方法，首先求解水流连续方程和动量方程，然后求解水流挟沙能力、泥沙不平衡输沙方程和泥沙连续方程，具体求解过程如下。

1.非恒定流计算

采用Preissmann隐式差分格式求解水流连续方程和动量方程，Preissmann格式网格布置见图2-10，因变量和其导数的离散格式为

式中：θ为加权系数，0≤θ≤1。

利用式（2-67）～式（2-69）可得到式（2-62）、式（2-63）的差分形式，对差分方程进行线性化，在线性化过程中，忽略增量的乘积项，最后得到以下线性方程组：

式中

以上式中：上角标为时间序列，下角标为断面序号；A_ij、B_ij、C_ij、D_ij、E_ij（i=1，2）为第j单元河段差分方程的系数（j=1，2，3，…，N-1，其中N为断面个数）。

式（2-70）、式（2-71）可针对任何一计算点（j，j+1）写出，如在模型中有N个计算点，就能对2N个未知数写出2（N-1）个这样的方程，再加上、下游两个边界条件，即构成2N个方程式组成的包含2N个未知数的方程组，因此方程组可以求解。

由于差分方程中的系数包含有未知数，方程求解不能直接求出未知变量，因此方程求解时必须进行迭代处理。下面给出求解的方法。

给定边界条件

对于方程组假定对于j点有如下线性关系：

可以证明在下一个点j+1也存在着类似的线性关系：

将式（2-73）代入式（2-70）、式（2-71）可得

其中

消去式（2-70）、式（2-71）两个方程中的ΔZ_j、ΔQ_j，然后将ΔQ_j+1表达为ΔZ_j+1的函数，得

比较式（2-74）、式（2-79）得

根据循环计算式（2-76）～式（2-78）和式（2-80）、式（2-81），在追的过程中求得系数H_j、I_j、J_j、F_j、G_j，然后根据式（2-73）、式（2-75），在赶的过程中求出ΔQ_j和ΔZ_j，进而计算出各河段的水位及流量。

2.支流汇入断面

对有支流汇入处，若汇流区附近支流水位与干流水位差别较小，则在汇入口上下各设一计算断面（图2-11），断面间距很小，为此可得如下关系式：

上式改写成

由于两个计算断面断面间距很小，则，上式进一步改写为

将ΔQ_j=F_jΔZ_j+G_j代入式（2-84）得

将式（2-86）与式（2-74）对比得

对照式（2-85）与式（2-75）可得

H_j=0，I_j=1，J_j=0

有了汇流口上下断面的H、I、J、F、G等参数，即可同其他断面一起应用追赶法计算各河段水力要素值。

3.不平衡输沙方程求解

利用迎风格式，将式（2-65）离散为差分方程，整理后得

当Q≥0时，利用上边界条件，自上而下计算各断面含沙量；当Q<0时，利用下边界条件由下至上计算各断面含沙量。

4.悬移质计算辅助方程

悬移质计算辅助方程包含挟沙力、挟沙能力级配方程及悬移质级配变化方程，在同样水流和床沙条件下输沙达到平衡时的悬移质级配，其方程可按三种输沙状态确定等均与韩其为的M1-NENUS-3模型相同，此处不再赘述。

（三）模型验证

采用2003年5月—2007年12月三峡水库实测水文泥沙资料对模型进行了验证。模型计算范围包括长江：朱沱至三斗坪；嘉陵江：北碚至重庆；乌江：武隆至涪陵。

五、一维非恒定水流泥沙数学模型（黄仁勇）

长江科学院黄仁勇提出的一维非恒定水流泥沙数学模型建立于三峡工程试验性蓄水阶段。应用于“十一五”三峡工程泥沙问题研究中的水库淤积计算^［18］。

（一）模型基本方程

水流连续方程

水流运动方程

泥沙连续方程

河床变形方程

以上式中：ω为泥沙沉速，角标i为断面号；Q为流量；A为过水面积；t为时间；x为沿流程坐标；Z为水位；K为断面流量模数；S为含沙量；S_*为水流挟沙力；ρ′为淤积物干容重；B为断面宽度；g为重力加速度；α为恢复饱和系数；A_d为河床冲淤面积。

（二）汊点连接方程

1.流量衔接条件

进出每一汊点的流量必须与该汊点内实际水量的增减率相平衡，即

式中：Ω为汊点的蓄水量，如将该点概化为一个几何点，则Ω=0。

2.动力衔接条件

如果汊点可以概化为一个几何点，出入各个汊道的水流平缓，不存在水位突变的情况，则各汊道断面的水位应相等，即

3.边界条件

计算中不对某单一河道单独给出边界条件，而是将纳入计算范围的三峡水库干支流河道作为一个整体给出边界条件，各干支流进口给出流量和含沙量过程，模型出口给出水位过程、流量过程或水位流量关系。

（三）模型求解

1.水流方程求解

采用三级解法对水流方程进行求解，首先对水流方程式（2-89）和式（2-90）采用普列斯曼的四点隐式差分格式进行离散，可得差分方程如下：

式中系数均按实际条件推导得出。

假设某河段中有mL个断面，将该河段中通过差分得到的微段方程式（2-95）和式（2-96）依次进行自相消元，再通过递推关系式将未知数集中到汊点处，即可得到该河段首尾断面的水位流量关系：

式中系数α₁、β₁、δ₁、θ_mL、η_mL、γ_mL由递推公式求解得出。

将边界条件和各河段首尾断面的水位流量关系代入汊点连接方程，就可以建立起以三峡水库干支流河道各汊点水位为未知量的代数方程组，求解此方程组得各汊点水位，逐步回代可得到河段端点流量以及各河段内部的水位和流量。

2.泥沙方程求解

对泥沙连续方程式（2-91）用显格式离散得^［42］

将式（2-91）代入式（2-92），然后对河床变形方程式（2-92）进行离散，得

式中：Δx为空间步长；Δt为时间步长；ΔA_di为河床变形面积；角标n为时间层。

在求出干支流河道所有断面的水位与流量后，即可根据式（2-99）自上而下依次推求各断面的含沙量，汊点分沙计算采用分沙比等于分流比的模式，最后根据式（2-100）进行河床变形计算。

（四）有关问题的处理

1.床沙交换及级配调整

关于床沙交换及级配调整，该模型采用三层模式，即把河床淤积物概化为表、中、底三层，表层为泥沙的交换层，中层为过渡层，底层为泥沙冲刷极限层。规定在每一计算时段内，各层间的界面都固定不变，泥沙交换限制在表层内进行，中层和底层暂时不受影响。在时段末，根据床面的冲刷或淤积向下或向上输送表层和中层级配，但这两层的厚度不变，而底层厚度随冲淤厚度的变化而变化。

2.水流挟沙力计算

水流挟沙力公式为^［43］

上二式中：p_L为第L组泥沙的级配；ω_L为第L组泥沙的沉速；S_*为水流总挟沙力；k为挟沙力系数，水库为0.03，天然河道为0.02。

3.恢复饱和系数α

恢复饱和系数是泥沙数学模型计算的重要参数，是一个综合系数，需要由实测资料反求。但是影响因素很多，既与水流条件有关，又与泥沙条件有关，随时随地都在变化，在大多数泥沙冲淤计算中都假定为一正的常数，通过验证资料逐步调整。该模型对泥沙冲淤采用分粒径组算法，如果对各粒径组都取同样的α值，由于各组间的沉速相差可达几倍甚至几百倍，因而从计算结果看，在同一断面上小粒径组相对于大粒径组来说其冲淤量常常可忽略不计，这往往与实际不尽相符。从三峡水库蓄水运用以来进出库的各粒径组泥沙实测资料来看，各粒径组泥沙的沿程分选现象均非常突出。目前有关恢复饱和系数α取值的研究有如下共识：①不同粒径组泥沙的恢复饱和系数值不同；②恢复饱和系数取值应随泥沙粒径的增大而减小；③恢复饱和系数值应随空间和时间而变化。为此本模型在计算中对不同粒径组泥沙的恢复饱和系数α_L采用如下经验公式计算^［44］：

上二式中：ω为干流计算进口朱沱站处的泥沙平均沉速，由于20世纪60年代系列和90年代系列进口朱沱站ω均与第5组泥沙沉速非常接近，为便于计算，最后统一取ω为第5组泥沙沉速ω₅；指数λ=c/J，其中J为水力坡度，c=0.833×10^-10，为坝址处多年平均流量。

4.糙率系数n的确定

糙率系数是反映水流条件与河床形态的综合系数，其影响主要与河岸、主槽、滩地、泥沙粒径、沙波以及人工建筑物等相关。阻力问题通过糙率反映出来，河道发生冲淤变形时，床沙级配和糙率都会作出相应的调整。当河道发生冲刷时，河床粗化，糙率增大；反之，河道发生淤积，河床细化，糙率减小。长系列年计算中需要考虑在初始糙率的基础上对n值进行修正。该模型根据实测水位流量资料进行初始糙率率定，各河段分若干个流量级逐级试糙。

5.节点分沙

进出节点各河段的泥沙分配，主要由各河段临近节点断面的边界条件决定，并受上游来沙条件的影响。该模型采用分沙比等于分流比的模式：

式中：Q_i，in、S_i，in分别为节点上断面流量和含沙量；S_j，out为节点下断面含沙量。

6.区间流量

三峡水库库区长度超过700km，区间支流众多，除嘉陵江和乌江来流外仍然存在着较大的区间流量，以干流朱沱站、嘉陵江北碚站和乌江武隆站三站进口实测资料统计，三峡库区的区间流量占出口宜昌站流量的百分比：20世纪60年代系列为12.2%；90年代系列为13.3%；2003—2007年为出口黄陵庙站流量的8.6%。区间流量往往会集中汇入，其影响不容忽视，在非恒定流计算中必须考虑区间流量的汇入。该模型将区间流量通过分配到各入汇支流上加入计算河段，各入汇支流流量根据进出库控制站及库区水文站已有实测水文资料通过计算得到。计算中没有考虑区间入库沙量。

（五）模型验证

采用三峡水库2003—2007年泥沙淤积实测资料和汉江丹江口水库汉江库区1968—1985年泥沙淤积实测资料进行了验证。

六、平面二维及一维嵌套泥沙数学模型（李义天）

在“七五”国家重点科技攻关项目“泥沙数学模型及糙率研究”子课题中，武汉水利电力大学李义天等建立了一维与二维嵌套泥沙数学模型，并在三峡工程变动回水区航道冲淤变化等泥沙问题研究中得到了应用^［20］。

（一）基本方程

由于所研究的问题是建立一维与二维嵌套数学模型，因此基本方程应包括一维模型及二维模型基本方程两部分。

1.一维模型

水流连续方程

水流运动方程

悬移质扩散方程

悬移质河床变形方程

推移质河床变形方程

水流挟沙力公式

推移质输沙率公式

以上式中：A为过水面积；Q为流量；v、H分别为断面平均流速和水深；Z为水位；n为糙率系数；R为水力半径；g为重力加速度；S、S^*分别为断面平均含沙量及水流挟沙力；N为悬移质粒径组数；、分别为悬移质和推移质淤积物干容重；Q_b为推移质输沙率；d为床沙粒径；η_sn、η_b分别为悬移质和推移质冲淤厚度，河床总冲淤厚度为。

2.二维模型

水流连续方程

水流运动方程

悬移质扩散方程

悬移质河床变形方程

推移质河床变形方程

水流挟沙力公式

推移质输沙率公式

以上式中：v、h分别为垂线平均流速和垂线水深；v_x、v_y分别为x和y方向的垂线平均流速分量；S_n、分别为垂线平均含沙量及垂线挟沙力；ε为紊动黏性系数；k为泥沙扩散系数；q_bx、q_by分别为x和y方向单宽推移质输沙率；为单宽推移质输沙率。

由于研究的问题是长时期的河床冲淤变化情况，因此在计算中略去了非恒定项［即∂A/∂t、∂（AS_n）/∂t、∂h/∂t、∂u/∂t、∂v/∂t及∂（hS_n）/∂t分别等于0］，亦即采用恒定流模型。

（二）定解条件

对一维模型，上游进口断面给定流量、悬移质含沙量及级配、推移质输沙率，下游出口给定水位，并且在河段内给定初始时刻河床地形及床沙级配。

对二维模型，在进口断面给定单宽流量、悬移质含沙量及级配、推移质输沙率沿河宽分布，在出口断面给定水位沿河宽分布，在河道两岸取v_x=0，v_y=0及∂S_n/∂n=0，∂q_b/∂n=0（∂/∂n为河岸边界的法向导数），并且在河段内给定初始时刻河床地形及床沙级配。

一维与二维嵌套模型，在一维模型和二维模型的连接断面上，按各模型对边界条件的要求，由一维与二维嵌套模型的连接条件确定。

（三）一维与二维嵌套模型的连接条件

在一维与二维嵌套模型的连接断面上，一维模型仅能给出断面平均量，而二维模型可给出沿河宽的分布量，一维模型和二维模型对边界条件的要求不同，因此存在着一维模型和二维模型的连接问题。按照连接断面所处的位置可分为两种连接断面：一种是上游为一维模型计算河段，下游为二维模型计算河段；另一种则是上游为二维模型计算河段，下游为一维模型计算河段。前者记为第Ⅰ种连接断面，后者记为第Ⅱ种连接断面。

1.水流运动连接条件

在恒定流模型中，水流运动的计算是从下游向上游逐段进行的。因此在连接断面上，由下游河段向上游河段提供水流参数的连接条件。

（1）水位。在连接断面上，水位应满足一般条件，即

式中：Z为断面平均水位；z为垂线水位；y为沿河宽方向坐标；B为河宽。

对第Ⅰ种连接断面，由下游河段的二维水力计算知道连接面上水位沿河宽分布，需确定连接断面上断面平均水位作为上游一维计算河段出口水位，由式（2-121）确定。

对第Ⅱ种连接断面，由下游河段的一维水力计算知道连接断面上平均水位，需确定连接断面上各垂线水位作为上游二维河段出口断面水位。如果连接断面选取在比较顺直、断面沿程变化较小的河段上，水位沿河宽变化不大，可取连接断面平均水位作为各垂线水位。

（2）流量及流速分布。在连接断面上应满足一般条件，即

流速沿河宽分布可以由水流平面图法给出，也可按由水流平面图法导出的垂线平均流速公式

给出，其中f（η）为确定糙率沿河宽分布的经验关系。

在第Ⅰ种连接断面上，由于水深未知，由式（2-123）可确定v、h关系，该式起边界条件作用，对第Ⅱ种连接断面，下游一维计算河段流量Q取模型进口断面给定的流量，由于在上游二维计算中流量始终也是取模型进口断面的流量，因此由上游二维计算得到的连接断面水深及流速分布可以满足式（2-122）。

2.泥沙运动连接条件

恒定流模型泥沙运动的计算是从上游向下游逐段进行的。因此在连接断面上，由上游河段向下游河段提供泥沙参数条件。

（1）悬移质含沙量连接条件。在连接断面上，悬移质含沙量应满足一般条件，即

对第Ⅰ种连接断面，由上游的一维泥沙计算知道连接断面平均含沙量，需确定连接断面上各垂线平均含沙量作为下游二维计算河段进口断面悬移质含沙量。根据水流挟沙特性，作为一种近似处理，假定各垂线平均含沙量与v³/（ghω）成正比，即

将上式代入式（2-124）确定系数C，然后将C代入式（2-125）得

对第Ⅱ种连接断面，由上游河段的二维泥沙计算知道连接断面上含沙量沿河宽分布，需确定连接断面平均含沙量，由式（2-124）计算。

（2）推移质输沙率连接条件^［45］。在连接断面上，推移质输沙率应满足条件

（四）有关公式的确定

1.阻力公式

一维过渡期糙率按式（2-128）计算：

式中：n为过渡期糙率；n₀为淤积前的天然糙率；n_k为淤积后最终糙率；a_k为与n_k相应的淤积面积；a为过渡期淤积面积。

二维糙率沿河宽分布按下式计算：

式中：n为二维糙率；n₀为一维糙率；J为二维比降；J₀为一维比降。

2.水流挟沙力公式

一维水流挟沙力按张瑞瑾公式计算^［41］：

分析实测资料得出如下二维床沙质水流挟沙力公式：

式中：K为断面平均水流挟沙力系数，其取值与一维水流挟沙力系数相同。

分组水流挟沙力的确定同时考虑了水流条件及床沙组成条件，对来沙条件的影响，则可通过冲淤造成的床沙级配变化来加以考虑。从这一认识出发，建立了如下床沙质水流挟沙力级配和床沙级配之间的关系：

式中：为第n粒径组泥沙的水流挟沙力级配；P_bn为第n组泥沙在床沙中所占比例；v_*为摩阻流速。

式中：σ_v为水流垂向紊速均方差；v′为水流的垂向紊动速度。

求出后，用乘以总水流挟沙力，即可求得分组水流挟沙力。式（2-132）曾用荆江资料进行了验证，结果符合良好。

式（2-132）用于天然情况下床沙组成属于卵石河床，修建水库后转化为沙质河床的冲淤计算时，需要确定床沙质的最小粒径（床沙质与冲泻质的分界粒径）及床沙质的最大粒径。

将阻力公式代入悬浮指标的表达式，并用实测资料确定有关参数，得到以下床沙质最小粒径及最大粒径计算公式：

与床沙质最小粒径对应的沉速

与床沙质最大粒径对应的沉速

3.推移质输沙率公式

在模型中选用了窦国仁公式^［41］：

式中：p_b为推移质粒径部分泥沙在床沙中占的比例；v₀为起动流速；ω为推移质平均粒径对应的沉速；k₁为系数，对川江卵石推移质k₁=0.01。

由于目前推移质实测资料较少，且精度较差，在一维及二维模型中采用了相同的计算公式。

（五）数值解法

数值解法与前述基本方程相对应，也分一维模型和二维模型数值解法两部分。在计算中为了节省计算工作量，均采用非耦合解，亦即先计算水力条件，再计算悬移质含沙量、推移质输沙率及河床冲淤变化等。控制时间步长，使在同一计算时段内，河床冲淤对水力条件的影响可忽略不计的假定得以成立。

1.一维模型

一维模型的水力计算相对比较简单，该模型采用了常用的推水面线方法。

含沙量采用下式计算：

式中：S、S^*分别为断面平均含沙量及挟沙力；Δx为上、下断面间距；q为单宽流量；α为系数；脚标带“0”表示为上游断面值，其余为下游断面值。

悬移质河床变形用以下差分方程计算

推移质河床变形用以下差分方程计算：

上二式中：Q_b0为上游断面推移质输沙率；Q_b为下游断面输沙率；Δt为时间步长。

计算床沙级配的方程为

式中：Δη_n为第n组泥沙冲淤厚度，包括悬移质和推移质；Δη为总的冲淤厚度；P_0bn、P_bn分别为时段初及时段末的床沙级配；h₀为床沙可动层厚度，其大小与河床冲淤状态、冲淤强度及冲淤历时等有关。

2.二维模型

（1）计算方法。二维泥沙数学模型计算方法用单元插值函数近似地代替其解析解，这样网格可划分成任意形状，与不规则岸线能较好地吻合。

单元插值函数采用九结点四边形单元（图2-12）进行插值，其形函数为

任意函数f可写成

式中：i为单元结点号；f_i为单元结点上函数值。

任意函数f的一阶、二阶导数为

其中

将f换成水深、流速等，即可求出对应的导数。

（2）水流计算。基本方程的离散是在计算区域内，对任一节点，将式（2-142）、式（2-143）代入式（2-113）、式（2-114）、式（2-115）得

以上式中：v_xm、v_ym、h_m分别为m节点流速、水深。

流场计算为先用一维方法推求水面线，用所得结果确定各断面水位的迭代初值。大致估算回流的范围。在主流区内用水流图法（累积流量法）估算各断面流速分布，并将其分解为x、y方向分速，用以确定主流区流速迭代初值。回流区流速迭代初值用下式确定：

式中：v_xb为靠近回流交界的流速在x方向的流速；a=1.5πB₁/B₂，B₁为从回流交界面算起的起点距，B₂为回流宽度。

用同样的方法可确定方向流速的迭代初值。

计算中采用了可同时消除数值波动的松弛迭代法：

式中：f为m点本次计算值；f_i为i点的前次迭代值；i为与m点相邻的结点号；N为与m点相邻的结点数；β为松弛因子。

对上述数值方法的精度及迭代顺序的收敛性等进行了研究，结果表明，在矩形网格的情况下，其精度与有限解析法、有限元法及有限差分法差别不大。在任意形状单元的情况下，窄长形单元误差较大，在划分单元时应引起注意。迭代序列在一般情况下是收敛的，在回流与主流交界面附近，有时会出现一定的问题，这一问题可通过适当选取松弛因子β加以解决。

（3）悬移质含沙量计算。在水流结构计算中，用水流平面图法求得流束分布迭代初值后，在主流区的每一条流束内均可用一维方法计算其相应的含沙量。其计算方法与前述一维模型中的式（2-136）相同，但式中符号的意义均为流速平均值，S^*用二维水流挟沙力公式式（2-131）计算。回流区含沙量的迭代初值，取含沙量沿河宽不变，其大小取靠近回流交界面流速的含沙量。

将式（2-142）、式（2-143）及相应的含沙量初值代入式（2-116），得在计算中，泥沙扩散系数k取值对回流区的淤积量有较大影响。经过检验取k=0.6v_*h。

（4）河床变形计算。求出各粒径组的含沙量后，用式（2-117）计算悬移质河床变形。将式（2-117）写成差分形式

将推移质输沙率公式式（2-135）改写成x，y方向的单宽推移质输沙率

式中：v为垂线平均流速。

将式（2-151）、式（2-152）及式（2-142）、式（2-143）代入推移质河床变形方程，并将时间导数项写成差分形式，得

河床总的冲淤厚度。床沙级配计算方法与式（2-139）相同。

（六）一维与二维嵌套模型的验证

利用1985年、1986年长江上游河段的水流和河道冲淤实测资料对模型进行了验证。