黄河下游河道准二维泥沙数学模型研究

摘要：本文首先简要回顾了以往黄河河道泥沙数学模型研究成果，再以水流运动方程及经过作者修正的泥沙运动方程为基础，同时引入与实测资料相符合的水流挟沙力、动床阻力、泥沙级配等计算公式作为补充方程，构造出黄河下游河道准二维泥沙数学模型。然后，采用1986年11月～1996年10月这10年长系列实测资料，开展了验证计算。其结果表明，该模型不仅能计算黄河下游河道一般洪水引起的河床冲淤变形，还能成功地模拟出大沙年下游处于强烈淤积时的规律。

关键词：黄河下游河道 准二维泥沙 数学模型

1 黄河河道数学模型研究的简要回顾

黄河数学模型的研究，是与流域规划、工程建设和管理运用等生产紧密结合的。早在1955年编制黄河综合利用规划技术经济报告中，就曾在黄河三门峡水库规划阶段用初级的一维恒定平衡输沙模型对水库淤积和下游河道的河床冲刷变形进行计算[1]。当时的计算结果认为，在桃花峪下游冲刷9年后河床刷深27m，冲刷量及冲刷速度显然比实际夸大甚多。

三门峡水库投入运用后，库区严重淤积，并迅速向上游延伸，与原来的计算结果截然不同。为研究改建方案，通过实测资料建立了许多泥沙冲淤量与水文要素之间的经验关系式。麦乔威、赵业安、潘贤娣等学者在分析下游河道冲淤变化及挟沙能力变化规律的基础上，提出了三门峡水库下游河床冲淤计算的方法[1]。

在1975年～1977年进行治黄规划过程中，为研究各种规划方案对黄河下游的减淤作用，黄河水利委员会及其有关科研、设计部门的科技人员，使用两种河道冲淤计算的经验模型开展了大量的计算工作，麦乔威、李保如两位学者进行了总结。其中，在规划方案的初步比较阶段，使用的是李保如提出的框算模型，只考虑影响输沙能力的最主要因素（来水、来沙）。对于问题比较复杂和初选的方案，则使用较详细的经验模型[2]。

20世纪80年代，刘月兰、韩少发、吴知等学者在上述经验模型的基础上，根据不同河段汛期 、非汛期输沙经验公式，引入水流连续方程、动量方程及泥沙平衡方程式，建立了一套颇具特色的黄河下游河道冲淤计算模型[3]。该模型在黄河规划中得到应用，只是在第三次加固大堤设计中，预计的下游淤积量与实际相差甚大。1987年～1990年，谢鉴衡、韦直林、魏良琰同刘月兰、张启卫等，开展了黄河下游洪水演进预报水动力学数学模型研究。特别是1992年之后，对国家“八五”重点科技攻关项目研究的全面启动，使黄河数学模型研究有了进一步发展的机遇，在钱意颖、张启卫主持下，开展了“黄河泥沙冲淤数学模型的应用”专题的研究[4]，在此期间，黄科院、清华大学、武汉水利电力大学及中国水科院的有关学者，对黄河数学模型的研究都有较大贡献。附带指出，黄委会有关专家曾对美国最通用的泥沙模型HEC－6模型，用黄河资料进行了检验，认为这类模型直接用于黄河必须做很多的工作。另外，美国学者杨志达博士等1986年开发了GSTARS模型，参照我国经验进行修正后，对黄河杨集至孙口河段进行了验证计算。以谢鉴衡为首的专家组的验收意见指出：“本模型要应用于包括游荡性河段在内的黄河下游，还需进一步做更多的研究工作”。此后，世界银行又委托沈学文对该成果进行评价，提交的书面意见认为：“模型可以用于黄河某些不平衡输沙问题不是十分严重的顺直河段，但横向输沙很麻烦且需要丰富的经验”显然，海外数学模型同黄河实际尚存在着很大的距离。

在黄河平面二维泥沙数学模型的研究方面，韦直林[5]采用有限元法对黄河下游夹河滩至高村河段的1982年的洪水演进过程进行了模拟：程晓陶、刘树坤[6]以二维恒定流与不平衡输沙理论为基础，将显式有限差分法与有限体积法相结合，并引入张红武水流挟沙力公式计算泥沙冲淤变化，研制出夹河滩至高村河段滩区水沙运行数学模型。结合工程实际，张世奇、杨国录等也对黄河京广铁桥以上河段的冲淤规律进行了数值计算。与此同时，张红武、江恩惠、赵连军、李炜、郑邦民、杨明、钟德钰、张凌武等学者联合攻关，开展了黄河下游不同河段的平面二维数学模型研究。例如，张红武、李东风、许雨新[7]等应用黄河泥沙运动基本理论研究的最新成果，采用有限元法建立了济南河段平面二维水流泥沙数学模型。从河床冲淤形态、流量水位过程、全域流速场等方面验证结果可以看出，该数学模型计算结果与实测资料符合较好。

应该指出，由于黄河问题的复杂性；现有的数学模型对所存在的参数，多进行经验处理或用实测资料率定。而黄河水沙条件及河床边界条件如果发生很大变化，那些经验性较强的参数的处理，往往会影响新条件下的计算结果，因而不能用来预测黄河未来冲淤趋势。另一方面，水流挟沙力及河床糙率是黄河泥沙数学模型应考虑的两个重要问题，不少模型中没有解决这两个问题。因此，即使在验证时差别不大，但在重要的特殊洪水期的计算结果，与原型却差异甚多，大大影响了数学模型在黄河防汛规划中的使用价值。

再者，恢复饱和系数的确定也是黄河泥沙数学模型的关键和难点所在，理论值应大于1，而一般的黄河模型中的取值都远小于1，其处理值得进一步商榷。为此谢鉴衡指出：“现阶段 对黄河数学模型上的恢复饱和系数，无论理论认识或实际经验都甚感不足，有必要在这两方面做进一步的工作。”我们认为，该问题主要是在假定非平衡状态下的恢复饱和

系数与平衡状态下相等之后出现的，存在着理论上的缺陷，因此，无论对该系数如何进行改进，都难以从根本上提高泥沙数学模型的理论水平。

作者针对上述存在的问题，长期开展了黄河数学模型的研制工作，本文介绍黄河下游准二维泥沙数学模型的研究结果。

2 准二维泥沙数学模型的建立

2.1 基本方程 本模型计算除选用描述河道水流运动的基本方程外，还改正了在泥沙数学模型中占据重要位置的河床变形方程及泥沙连续方程式[8，9]，即：

河床变形方程

泥沙连续方程

式(1)～(2)中：i为断面号；j为子断面号，河床高程最低的子断面号j为1，最高的取j为m；A为过水面积；t为时间；x为沿流程坐标；Z为水位；ω为泥沙浑水沉速；S为水流含沙量；S为水流挟沙力；γ0为河床淤积物干容重；bij为子断面宽度；Zbij为子断面平均河床高程；α为平衡含沙量分布系数，由下式计算[10]：

其中：

式中：κ为卡门常数；cn为涡团参数(cn=0.375κ)；C为谢才常数；u为摩阻流速；且有：

运用量纲分析及相似转化原理，并根据动床河工模型试验结果，可求得Ks的表达式为：

对于非饱和系数fs，通过归纳分析，再借助试验资料，得出：

2.2 水流挟沙力及糙率计算 为保证数学模型适用于黄河下游高含沙洪水，提高多沙河流输沙计算的精度，水流挟沙力计算采用黄河上常用的张红武水流挟沙力公式[10]

式中：

式中：ωsk为第k组粒径泥沙在浑水中群体沉速；ωok为第k组粒径沉速；D50为床沙中值粒径；d50为悬沙中值粒径(mm)；psk为第k组粒径泥沙的重量百分数；N为非均匀沙的分组数。

经舒安平、江恩惠、朱太顺及陈雪峰、陈立、李义天[11]等学者，通过黄河、长江等河流大量实测资料系统检验结果表明，式(8)是现有公式中最适用于黄河等天然河流的水流挟沙力公式，目前已得到广泛应用。而且教科书中也被推荐为适用于高含沙水流的公式[12]。

本模型采用作者的动床糙率公式开展黄河洪水水位模拟计算，可较好地描述水力泥沙因子的变化对摩阻特性的影响，并能反映天然河道中各种附加糙率的影响。该糙率公式为：[13]

式中：δ为摩阻厚度，采用下游资料，得出δ与床沙中径D50及佛汝德数Fr的关系式为：

2.3 悬移质泥沙与床沙级配的计算 从悬浮于水体中泥沙颗粒的受力分析入手，经理论推导和经验处理可建立悬移质泥沙级配及粒径二阶圆心距的计算公式[14]

上述3式中，P(di)为小于某粒径的泥沙所占的百分数；d50为悬沙中值粒径；dcp为悬沙平均粒径；m为指数；ξd为悬沙粒径的二阶圆心距，其数值可表征出泥沙组成的非均匀程度。由式(15)、(16)可知，对于变量d50、dcp、ξd，已知任意两个，即可确定第三个变量的值，悬沙组成即可由式(14)求解出来。同时，通过黄河下游 水文站实测资料回归分析，建立了悬沙中值粒径d50和平均粒径dcp之间的如下关系式：

式中：vm为浑水运动粘滞系数。

对于天然细沙河流，底沙随时可能跳起，同悬沙的交换几率很大。我们以跃动模式分析床沙起动时的受力入手，推导出与悬沙相类似的级配及粒径二阶圆心距公式：

式中：P(Di)为床沙中小于某粒径的泥沙所占的百分数；Dcp为床沙平均粒径；m为指数；ξD为床沙粒径的二阶圆心距。

从悬沙与床沙交换过程中任一粒径组的沙量平衡入手，经理论推导建立了如下悬沙与床沙平均粒径与二阶圆心距变化规律的表达式：

式(21)～式(24)中：dL、ξL分别为侧向入流泥沙的平均粒径及粒径的二阶圆心距，P为床沙混合层厚度，dc、ξc为冲淤物平均粒径与粒径的二阶圆心距，dc与P的具体计算方法

通过对黄河下游实测资料回归得出如下悬沙平均粒径沿横向分布公式：

式中：dcpi为断面平均悬沙平均粒径，dcpij为断面上任一点悬沙平均粒径。C2为断面形态系数，由沙量守衡可求得

2.4 悬移质含沙量横向分布计算 通过对黄河大量实测资料分析发现，含沙量的横向分布规律不仅与水力因子、含沙量大小有关，还与悬沙组成密切相关。悬沙粒径越细，含沙量的横向分布越均匀，为此除引入含沙量因子外，还引入悬浮指标ω/ku*，建立了含沙量横向分布公式，即

式中：Vi、Vij分别为断面平均及任意一点的流速；Hi、Hij分别为断面平均及任意点的水深；u为断面平均摩阻流速；C1为1左右的断面形态系数，由沙量守恒可求得：

式中：qij为断面任一点单宽流量，y为横向坐标；a、b为断面河宽两端点起点距。

2.5 河槽在冲淤过程中河宽变化规律的模拟 本数学模型所采用的模拟河槽宽度变化过程的方法，是根据计算河段应受特定的河相关系均衡调整原理进行的。模拟河段内各横断面宽度修正值，除受到给定河岸条件的限制外(例如受到抗冲河岸、山嘴及河道整治工程等的限制)，还要随着河流造床过程而自动调整。本模型选用张红武河床综合稳定性指标[9]作为河相关系均衡调整准则计算，即：

式中：B为河宽；i为河床纵比降。作为初步探讨，将河床的均衡关系式运用于河床断面的自动调整过程，对于河宽(包括子断面宽度)的调整以式(29)的约束关系确定其变化值。

3 验证计算结果

该数模建立之后，主要针对黄河下游诸如“82.8”、“77.7”、“92.8”、“96.8”等不同类型洪水开展了验证计算①。为进一步检验该数学模型对不同水沙系列年的适用性，特选用三门峡水库蓄水运用及滞洪运用初期1960年9月～1965年6月的水沙系列和近期的1986年11月～1996年10月水沙系列以及相应的冲淤变形资料，对模型进行检验。前者见文献[15]，本文只介绍后者的验证计算结果。

1986年～1996年，黄河下游总来水量为2979.6亿m3，总沙量为85.742亿t，其中汛期来水量为1404.8亿m3，占总来水量的47.1%，来沙量为81.7亿t，占总来沙量的95.3%。特别是1988年、1992年、1994年、1996年，这4年为丰沙年，汛期小浪底水文站平均含沙量均超过80kg/m3，汛期来沙量为48.6亿t，占总来沙量的56.6%。由于水沙搭配不合理，造成了下游河道严重淤积，10年内下游河道共淤积25.0亿t，其中铁谢至高村河段淤积19.2亿t，占总淤积量的76.8%，高村至艾山河段淤积3.6亿t，占总淤积量的14.4%，艾山至利津河段淤积2.2亿t，占总淤积量的8.8%。从淤积的年际年内分配上看，河道主要在汛期淤积，共淤积31.4亿t，而非汛期河道基本处于冲刷状态，共冲刷6.4亿t。其中1988年、1992年、1994年、1996年四年汛期淤积最为严重，共淤积22.4亿t，占总淤积量的89.5%。从泥沙的淤积部位看，大量泥沙主要淤积在主槽内，引起河槽高程明显抬升，平滩流量大幅度下降。如“92.8”、“96.8”两场洪水，花园口站洪峰流量分别为6020m3/s、7600m3/s，接近于该站洪峰的多年平均值，仅相当于中常洪水，花园口站却先后出现有水文记载以来的最高水位，洪水漫滩范围及传播时间也明显增大，洪水造成的经济损失巨大。

采用小浪底站实测水沙过程及悬沙组成作为进口边界条件。伊洛河、沁河入黄水沙及泥沙级配过程分别采用黑石关、小董(或武陟)站实测资料，沿黄引水引沙量采用实测资料沿程分配给出，作为已知的分流条件。初始地形采用黄河下游1986年汛后92个实测大断面控制，并将每一个大断面概化为12～15个子断面。初始床沙组成根据花园口、夹河滩、高村、孙口、艾山、泺口、利津等7站的实测床沙组成按断面间距内插或外延。出口边界水位采用相应年份 利津站水位流量关系控制。

表1给出了铁谢至利津不同河段不同时段计算和实测的冲淤量。计算表明，1986年11月～1996年10月，铁谢至高村河段淤积13.588亿m3，高村至艾山淤积2.445亿m3，艾山至利津淤积1.654亿m3，与相应的实测值13.707亿m3、2.581亿m3、1.571亿m3接近。

从上述计算结果可以看出，无论河道处于冲刷、淤积或冲淤交替状态，即使是遇到大沙年(包括其中的高含沙洪水)处于强烈淤积时，该模型均能很好地模拟黄河河道的输沙特性及河床变形规律。同时，该模型对河道在不同水沙组合条件下发生的上冲下淤、上淤下冲及沿程冲刷或淤 积等规律的模拟精度也令人满意。

4 结 语

本文首先回顾了前人关于黄河河道泥沙数学模型研究成果，再以黄河泥沙基础理论最新成果 为基础，构造出黄河下游河道准二维泥沙数学模型。然后采用原型实测的水沙过程及相应的河床变形资料，开展了验证计算。结果表明：该模型不仅能模拟下游河道长系列年一般洪水引起的河床冲淤变形，还能成功复演大沙年下游处于强烈淤积时的实际情况。该数学模型已在黄河下游防洪及规划中应用，而且还在小浪底水库运用方式研究中发挥重要作用。目前被黄 河主管部门确定最为可靠的黄河下游泥沙数学模型。

参 考 文 献：

[1] 麦乔威，赵业安，潘贤娣.多沙河流拦沙水库下游河床演变计算方法[J].黄河建设.1965，(3).

[2] 李保如，华正本，樊左英，陈上群.三门峡水库拦沙期下游河道的变化[c].见：第一次河流泥沙国际学术讨论会论文集(第一卷).北京：光华出版社.1980.

[3] 刘月兰，韩少发，吴知.黄河下游河道冲淤计算方法[J].泥沙研究.1987，(3).

[4] 钱意颖，张启卫，李松恒.黄河泥沙数学模型的研究与应用[C].黄河泥沙.郑州：黄河水利出版社.1996.

[5] 韦直林.黄河下游水流泥沙平面二维数学模型研究[J].武汉水利电力学院学报，1990，6：77-86.

[6] 刘树坤，宋玉山，程晓陶，等.黄河滩区及分滞洪区风险分析和减灾对策[M].郑州：黄河水利出版社，1999.

[7] Zhang Hongwu， Li Dongfeng， Xu Yuxin， Zhang Junhua. A Two Dismesional Sediment Mathematical Model for Unsteady Flow in the Lower Yellow River [J]. International Journal of Sediment Research.1999，14(2).

[8] 张红艺，杨明，等.高含沙水库泥沙数学模型的研究及应用[J].水利学 报，2001，(11).

[9] 张红武，吕昕.弯道水力学[M].北京：水利电力出版社，1993.

[10] 张红武，江恩惠，等.黄河高含沙洪水模型的相似律[M].郑州：河南科技出版社 ，1994.

[11] 陈雪峰，陈立，李义天.高中低浓度挟沙力公式的对比分析[J].武汉水利电力大学学报，1999，(5).

[12] 张瑞瑾，谢鉴衡，等.河流泥沙动力学[M].北京：水利电力出版社，1998.

[13] 赵连军，张红武.黄河下游河道水流摩阻特性的研究[J].人民黄河，1997，(9).

[14] 赵连军，张红武，江恩惠.冲积河流悬移质泥沙与床沙交换机理及计算方法研究[J].泥沙研究，1999，(4).

[15] 张红武，黄玉婕，赵连军，等.小浪底水库运用初期下游河道冲淤数学模型预测计算[J].水力发电学报.2002，(1).