摘 要 本文介绍的黄河下游宽浅河段汛期水沙运动仿真模型建立在二维不恒定流与不平衡输沙理论的基础上，采用无结构不规则网络，结合了有限差分法和有限体积法的优点，能较充分地反映滩区、生产堤与滩槽冲淤变化对洪水演进的显著影响，并具有实时动态图象显示和人机对话功能.
关键词 二维不恒定流，水沙运动仿真.
黄河下游滩宽槽浅，汛期来水来沙变幅很大.大堤间滩地面积占河道面积83%，有居民近140余万.滩区生产堤纵横，在纵比降为1/6000的河道内，滩面横比降可达1/2000；平滩流量下降到约4000m3/s.洪水一旦破堤漫滩；顶冲大堤，顺堤行洪，甚至倒灌，对滩区人民生命财产及黄河大堤的威胁更大.漫滩行洪和滩槽水沙激烈交换可引起洪峰严重变形，且河槽冲淤幅度与水位涨落达到同等量级.大幅快速的河床调整又反过来影响水流的挟沙能力与河道的行洪能力，引起水位-流量关系的复杂变化，使洪水预报异常困难.
对这种大范围的、在计算域内部受堤防溃决等突发性因素影响十分显著、水沙交换激烈的动边界、动床洪水现象，一般一维或一维扩展模型很难模拟.已有的洪水预报模型或河势演变模型对漫滩洪水大都作了简化，难以充分体现生产堤溃决、滩地行洪与滩槽水沙交换对洪水传播变形的影响.新研建的黄河下游河道水沙运动仿真模型，既要求计算精度高，以便为滩区洪水风险分析提供科学依据；又要求运算速度快，在大洪水发生时能及时准确为滩区居民发布洪水预警报信息.如此矛盾的要求增加了模型设计的难度.
目前，黄河下游夹河滩至高村(包括封丘倒灌区)防洪险要河段的水沙运动仿真模型已通过鉴定，其主要特点包括：①以一、二维不恒定流与不平衡输沙理论为基础，采用无结构不规则网格，可以较好地概化河道主槽、滩地的地形特征，充分反映生产堤、渠堤等对水沙运动的影响.②泥沙计算在黄金池方法的基础上作了改进，引入了张红武的一套挟沙能力公式［1］，使黄河下游主槽冲淤变化对洪水演进的影响能够较好体现出来，并能适应高含沙洪水的计算.③模型具有在运算过程中实时动态显示流量、水位、洪水淹没范围或泥沙冲淤计算成果的功能.④模型中设计了人机对话的窗口，对生产堤在洪水期间是否发生决口这种带有随机性的问题，模型使用者有机会根据实际情况或经验进行选择.⑤模型运算速度显著提高，对1427.8km2的计算域模拟6天洪水过程，仅需十余分钟机时，可满足为洪水预报服务的要求.⑥开发了模型的菜单对话式图形信息管理系统，使模型所需大量基本信息的前处理与计算成果的后处理工作大为简便.⑦计算结果与损失评估模型相结合，可以快速算出相应洪水条件下的受灾人口及可能造成的经济损失.
1 水沙运动的计算方法
1.1 基本控制方程 从黄河水沙运动特点出发，模型中选垂向平均单宽流量Q与水深H作为基本状态变量，其深度积分形式的控制方程组可记作如下的质量与动量守恒形式：
水流连续方程
(1)
水流动量方程
(2)
(3)
式中：H—水深；Z—水位；q—源汇项，本模型中代表有效降雨强度，不计降雨情况下q=0；M，N—分别为x、y方向的垂向平均单宽流量；u,v—分别为垂向平均流速在x、y方向的分量；n—曼宁糙率系数；g—重力加速度.
基本状态变量与其他交量之间有以下转换关系：
Z=H+B，
(4)
u=M/H，
(5)
v=N/h，
(6)
式(4)中B为地面高程.
在考虑泥沙冲淤计算时，基本控制方程组中要增加泥沙运动方程与河床变形方程.根据不平衡输沙理论，这两个方程通常可表达为如下形式：
泥沙运动方程
(7)
河床变形方程
(8)
式中：S,S*—分别为深度平均含沙量与水流挟沙力；ω—悬移质泥沙深度平均沉速；α—悬移质泥沙恢复饱和系数；es—悬移质泥沙扩散系数；γ′—淤积物干容重；其余符号与水流模型意义相同.
1.2 水流控制方程的简化与离散方式 黄河滩区计算域大，河势复杂多变，控导工程、生产堤、渠堤密布其间.在规则网格情况下，若为省计算时间而将网格空间步长放大，则会引起严重的地形失真，难以满足洪水风险评价对计算精度的要求.若加密网格以提高对地形的分辩率，又必然会增大计算量，难以满足缩短运算机时、为洪水预报服务的要求.
为此，本模型在开发中保留了体积积分形式的有限差分法的优点［2］，如在网格周边通道上计算流量、泥沙通量，在网格的形心计算水深、含沙量、挟沙力与底床的冲淤变化；同时又吸收了有限体积法中非结构性不规则网格的优点［3］.通道上定义的单宽流量，不限于正交的X或Y的方向，而可以是任意的通道法线方向的流量，同时，水深与流量在时间轴上分层布置，交替求解，物理意义清楚，并有利于提高计算稳定性，如图1所示.
a.H与Q的空间布置方式 |
b.H与Q的时间交错计算方式 |
图1 基本状态变量的离散化布置方式 将式(5)，(6)代入式(1)，则水流连续方程式改写后对计算域进行面积分，再利用高斯定得，可得
(9) 式中：u为计算域边缘上任一点的流速矢量，n为该点的外法线方向单位向量.令
(10)
对任一K边形网格，上式左端第2项的线积分可记为：
(11) 式中，A为网格面积，L为通道长度，下标k为对K边形网格逆时针方向的通道编号.
(12) 式中下标i为网格的编号.
(13) 式中hj为j通道上的平均水深，一般由下式计算：
(14) 式中DLj=DCj1+DCj2，为空间步长，如图2(a)所示.
(15) 式中Bm为相邻网格中较高的河底高程.
(16) 式中DL的定义为相临网格形心距在通道法线方向的投影，如图2(b)所示. |
a.河道型通道 |
b.滩地型通道 |
图2 不规则网格空间步长的分类处理 为适合黄河的特点，模型开发中对水流与泥沙计算的结合、网格形心流速的算法、滩唇的处理方法、出口断面无水位-流量关系及滩区横比降引起的数值效应等问题均作了有益的探讨，受篇幅限制，详细计算与处理方法将另文介绍，或请参见文献［4］. 2 模型的验证 |
图3(a) 1982年高村站计算与实测流量过程比较 |
图3(b) 1982年高村站计算与实测水位过程比较 |
| 布合理(图5)；④滩槽冲淤分布合理，充分体现了大洪水期间“淤滩刷槽”的现象，主槽中出现了一条清晰的连续的冲刷槽，反映了黄河汛期“大水出好河”的特点.根据实测资料，“82.8”洪水期间，泥沙淤积主要分布在生产堤间嫩滩上，约占滩地淤积量的95%，而生产堤外滩地的淤积量很少，计算结果表明泥沙淤积也是绝大多数分布在生产堤间的滩地上.主槽与嫩滩冲淤幅度的量级是合理的，主槽冲刷深度约为1.0—2.5m，嫩滩淤积厚度约为0.2—1.0m. |
图4(a) 假定生产堤不溃的流量过程模拟结果 |
图4(b) 假定生产堤不溃高村站水位过程的比较 |
图5(a) 1982年洪水黄河北岸大堤计算最高水面线与实测资料比较 |
图5(b) 1982年洪水黄河南岸大堤计算最高水面线与实测资料比较 |
3 模型的实时动态显示功能 以往的数值计算模型总是先“盲算”，只有待计算结束后，再进行数据分析，才能知道计算结果.而黄河水沙运行模型计算域广、数据量大，为了改变模型调试、运行效率低的状态，并适应计算中人-机对话的需要，专门开发了一整套绘图子程序，与计算程序一起运行，使各种计算结果能够实时动态显示出来，包括：①以不同颜色表示的水深或泥沙冲淤深度的平面分布.两者之间在计算中可以随时用指定的功能键进行往复切换，想看水深分布还是冲淤分布，随意选择.②以不同颜色线段表示的若干指定河道断面的流量变化过程与若干指定地点的水位变化过程.在有实测资料的情况下，还可以显示计算结果与实测过程的比较。③一些文字信息，如水量平衡系数、计算时间步长、计算时间步数、以年月日时表示的计算时间等.④平面图的比尺.当计算域调整时，只需改变图的中点坐标和比尺系数，即可调整模型的位置和大小，同时显示的平面图的比尺也自动调整. 4 模型的人机对话功能 5 结束语 参考文献 1 张红武，等.黄河高含沙洪水模型相似律的研究.黄委会水科院，1993，12. 2-Dsimulation model of flood routing and sediment transport in the Abstract Based on the theory of unsteady flow and non-equilibrium sediment transport, a 2-D simulation model of flood routing and sediment transport in the lower reach of the Yellow River has been developed. The model is built by non-structural irregular grids, combining with finite difference and finite volume methods, in which the significant influences of overbank storage, local leve and scouring (silting) in the process of flood propagation can be well considered. |
[责任编辑:simuhxf]
，则Q为任意n方向垂向平均单宽流量，当n取为x、y方向的单位向量时，Q即为M、N.
