关键词 黄河下游，水沙运动，数值模拟。

黄河下游花园口～孙口水沙运动仿真模型，是在“八五攻关”研建的夹河滩～高村模型的基础上，进一步向上下游扩展的成果[1，2]。模型计算范围包括了花园口至孙口河段黄河大堤间的主槽、滩地及封丘倒灌区，河道总长319km，计算域总面积达2340km2(图1)。

图1 黄河下游花园口～孙口河段二滩水沙运动仿真模型整体计算范

模型扩展面临两大难题。一是模型如何适应花园口～孙口河段河势特性的显著变化；二是模型如何适应黄河来水来沙变幅大的特点。新扩展的花园口～夹河滩河段为游荡性河段，主槽顺直，摆幅大；高村～孙口为游荡性向弯曲性过渡的河段，主槽弯窄，滩地小而零散。而原夹河滩～高村为有控制的游荡型河段，主槽摆幅较小，左右两岸滩地大而完整。为此模型扩展中从增强泥沙冲淤计算的自适应能力入手，采取了分段建模的技术路线，分段研究计算方法与参数选择的规律，最终形成了三个独立的分段模型与一个联接分段模型进行总体计算的整体模型。通过对“96.8”与“82.8”两场洪水的反复验证，不仅在模型的计算方法上有所改进，对有关参数的作用和组合规律也进行了初步的探讨，使模型向实用化目标又迈进了一大步。

1 模型总体结构设计

1.1 扩展模型的构造方式 在建成花园口～夹河滩、夹河滩～高村、高村～孙口三个分段模型之后，整体模型的开发考虑了“硬连接”与“软连结”两种构造方式。

“硬连接”保持运行模型结构不变，开发了专用的模型拼接软件，可将三个分段模型的基础信息数据集连接成一个整体进行计算。但是模型扩展后计算域又长又大，模型的动态图象演示与人 机对话功能受到影响，也不便于在模型中分段反映河势的差异。

“软连结”是保持三个分段模型的基础信息数据集不变，而对计算模型的结构作彻底改造，同样达到整体计算的目的。计算中模型可以通过功能健的切换，分段清晰地动态演示水沙运动过程，保留了通过人 机对话处理堤防决口或抢险堵口加高等现象的优点，并可以分段调整泥沙冲淤的计算方法与参数。计算结果亦分段保存，便于在防洪会商系统中利用已移植的图形软件分段演示计算结果。

图2 黄河下游花园口～孙口河段二维水沙运动仿真模型总体结构框图

1.2 扩展模型的结构框架 扩展模型开发期间面对面两大变化。一是外部运行环境的变化。黄委会的计算机系统开始联网，模型能够直接从网络上获取水情沙情信息作为计算边界条件，并需要将计算结果直接通过网络传送到防洪大厅，供决策会商系统调用；二是内部运行方式的变化。在采用“软连结”方案的情况下，模型系统内部的数据管理任务加大，同时增添了一些新的基础信息与计算信息的前后处理工作。扩展模型的结构框图如图2所示。

1.3 扩展模型的三种运行方式

由于黄河下游汛期河床变形剧烈，洪水的水位流量关系复杂，传统洪水预报方法有一定局限性，将水动力学模型用于黄河下游花园口至孙口河段的洪水预报，在80年代后期就曾有过积极的尝试[3]。本次将仿真模型扩展到花园口～孙口河段后，模型需要并有可能为洪水预报服务。历史洪水重现计算、防洪预案计算与洪水预报三种情况下，模型对水沙边界条件的要求差异较大(表1)，因此，运行方式也有所不同。模型三种运行方式的比较如图3所示。

图3 模型三种运行方式的比较

2 泥沙计算子模型的改进

2.1 泥沙计算子模型的特点与存在的问题 黄河下游二维水沙运动仿真模型中，泥沙运动子模型是在黄金池方法的基础上开发的[4]。模型中水流挟沙能力的计算，改用了张红武的一套公式[5]，较好反映了黄河汛期泥沙多来多排的特点。但是为了减少计算量，泥沙计算尽可能作了简化。模型上下延伸后，河道长达319km，上段河道泥沙计算失真对下段河道计算的影响加大。尤其对“96.8”洪水来说，含沙量变幅很大，迫切需要增强模型对泥沙计算的适应能力。

2.2 泥沙计算子模型改进的主要方面

(1)泥沙运动方程的求解中，进一步考虑了扩散项的影响。泥沙运动方程为：

（1）

对于不规则网格，式(1)中具有空间二阶导数特点的扩散项难以直接离散求解。因此，先仅对等式左端进行面积分：

（2）

式(2)经过简化离散后，可先求得网格含沙量的第一分步的中间计算结果Si1:

（3）

式(3)整理后，得

（4）

式中：Hi0—泥沙计算前一时段网格i的水深；Hi1—当前时段网格i的水深；DTs—泥沙计算的时间步长；Si0—网格内前一时段水流含沙量；Si1—网格含沙量第一分步的中间计算结果；Ai—网格I的面积；—DTs时段内通过Lik的平均单宽流量。Lik—i网格k通道的宽度；Sik—网格间通道上的含沙量；然后，由式(1)右端的扩散项求网格含沙量第二分步的中间变量Si2.令：

（5）

将式(5)右端改写成拉普拉斯算子的形式并作面积分，得：

（6）

对上式左端，假定HS的时间导数在网格中是均匀分布的；对上式右端，利用格林公式将面积分转换为线积分，得：

（7）

式(7)中n表示边线的外法线方向。该式的离散化形式为：

（8）

式中：Si2—i网格含沙量第二分步的中间计算结果；—i网格的平均水深，=(Hi0+Hi1)/2；DLk—k通道两侧网格的形心距；正、负号选择的规定—当下标i=k2时取‘+’；当i=k1时取‘-’；其余符号与上、下标意义同前。

(2)水流挟沙能力公式的实用形式。天然河道在不断变化的来水来沙条件下，通过河床冲淤变化，不断地调整自身的输水输沙能力。黄河下游河道作为多沙河流，洪水过程中这种短期的自调节现象更为显著剧烈[6]。在河流水沙运动仿真数学模型中，水流挟沙能力不仅是河道冲淤状态的重要的判别条件，而且是左右冲淤速率和幅度的重要因素。

从实用的目的出发，模型改造中，不再将影响因素考虑更细，而是设法增强模型的自调节能力。改造后的模型，对张红武挟沙能力公式中的泥沙沉速、粒径、级配仅采用一定时段的平均值，计算中允许在入口边界上分时段进行调整；同时将公式左端系数改造成冲淤量的函数，使挟沙能力随淤积而增强，随冲刷而减弱，形成负反馈机制。

同时，泥沙子模型的计算步骤进一步优化，改进后的运行框图如图4所示。

图4 泥沙子模型运算框图

3 模型人机对话功能的改进

八五攻关期间，模型的人机对话功能主要解决了堤防不同情况下的溃口问题。96.8洪水期间，滩区群众根据经验，对预报的7000～8000m3/S的洪水，采用了积极的抢险活动，生产堤被临时加高，包括一些预留的生产堤口门，也被堵了起来。为此，在模型的自动预警式与强制式对话中均增加了堵口与堤防加高的功能，同时增加了对话中误操作的识别功能，消除了以往由于误操作而导致计算意外中断的现象。新的自动预警式与强制式人机对话过程如图5、图6所示。

图5 自动报警式人机对话流程

图6 强制式人机对话流程

4 结语

　　在“八五攻关”成果的基础上，黄河下游二维不规则网格水沙运动仿真模型已扩展到了花园口～孙口河段，模型对泥沙计算方法有所改进，信息处理功能进一步增强。模型以96.8与82.8两场洪水作了验证，计算结果与实测资料吻合良好(见文献[2]).通过验证，对黄河下游花园口至孙口河段不同特性洪水的主槽演进、漫滩过程及冲淤特性有了更深入的了解。在验证中发现泥沙冲淤计算的合理性对于验证结果的影响十分显著。由于泥沙计算方法的改进，模型不仅对82.8洪水演进模拟效果良好，滩槽水沙交换与泥沙冲淤分布的计算结果更为合理，而且对96.8洪水在传播中的异常情况，尤其是对孙口断面“两峰合一”的现象也能合理地进行模拟。模型建成之后参加了1997年黄河防汛演习，在获取花园口预报流量过程之后，可迅速算出洪水漫滩演进过程。经与黄委中芬项目办合作，计算结果可以直接通过网络传输到防洪会商系统，在大屏幕上显示洪水泛滥过程及统计信息，供决策者参考。在98年黄河下游挖河固堤启动工程前期工作中，利用该模型对黑岗口～柳园口段的挖河情况进行了计算。计算结果与物理模型实验结果基本一致，并对挖河与修建控导工程对调整河势的作用、新河围堰破口过水时机、挖河宽度等提出了非常有价值的建议。

作者单位：1中国水利水电科学研究院 2黄河水利委员会河务局

致谢 宋玉山、胡一山、张红武、刘月兰、陈吉康、刘树坤等专家对模型的改造提出了宝贵的意见，仅此深表谢意。

本文于1998年11月11日收到。

参 考 文 献

1 程晓陶，陈喜军，薛云鹏等.黄河下游滩区水沙运行数学模型及减灾措施的研究，八五攻关子专题研究报告(85-926-01-02-01).1995，12.

2 程晓陶，蓝 虹，薛云鹏等.黄河下游水沙运动仿真模型扩展研究，中国水利水电科学研究院研究报告.1997，12.

3 韦直林等.黄河下游洪水预报水动力学数学模型，“黄河下游变动河床洪水演进与洪水位预报水动力学数学模型研究”报告之一，武汉水利电力学院河流数值模拟研究室.1990，3.

4 黄金池等.东平湖分滞洪区水沙运行数学模型及减灾措施的研究，八五攻关子专题研究报告.1995，12.5 张红武，江恩惠等.黄河高含沙洪水模型的相似律.河南科学技术出版社，1994，1.

6 钱宁，张仁等.黄河下游挟沙能力自动调整机理的初步探讨，地理学报.1981，36(2)