MacCormack法应用于输沙计算的研究

高波 黄金柏

摘 要：本研究探讨将MacCormack方法应用于河流的推移质计算，为流域输沙计算提供科学的数值方法。研究基于MacCormack法的差分方案，将水流计算分为预测计算和修正阶段，而后将推移质输沙计算公式与MacCormack计算方案耦合，并应用于设定条件的事例计算。计算结果形象地再现了河流在改变流量条件下输沙量随时间的变化，该方法适用于计算常流和变流量条件下的河道输沙量、也可用于河床演变，河道疏浚的计算。

关键词：MacCormack法，有限差分，推移质，输沙

中图分类号：TV148文献标识码： A 文章编号：

1前言

流域的水沙过程包括降雨产流、汇流的全过程, 以及泥沙侵蚀、输移、沉积的不同运动形态. 对这一物理过程的全面描述有其必要性和现实意义[1]。输水的同时输沙, 对河流起着泄洪排沙、维持河道正常演变的作用[2]。泥沙运动既应满足一定的力学规律, 同时也是一种随机现象[3]。输沙计算，一般先求解水流条件, 得到水深、流速、流量等水力因子, 将水力因子代入输沙方程中, 得到含沙量或输沙率过程[4]。

在现代河流泥沙物理模型及数学模型中,经常遇见输沙计算问题, 但由于河流动力学所面临的自然现象, 是边界复杂多变的天然河流及水流泥沙运动的多相体, 许多问题的运动规律尚无定论, 特别是水流挟沙力的计算, 人们还不能全面地把握它的力学机制, 因此, 多少年来, 泥沙界的学者们对水力输沙特性在不懈地探索[5]。崔侠等（1987）从冲淤过程的机理出发, 推导出了二度恒定非均匀流中泥沙扩敬方程, 求出了含沙量沿程变化规律, 着重讨论了水流挟沙力与含沙量的关系[6]。张洪武等（1992）从水流能量消耗和泥沙悬浮力之间的关系出发, 考虑了泥沙存在对卡门常数和泥沙沉速等的影响, 给出了半经验半理论的水流挟沙力公式并在黄河的输沙计算中应用[7]。胡海明等（1996）较为系统地研究了非均匀沙运动机理及输沙率的计算方法, 在建立泥沙运动交换模式的基础上, 指出现有非均匀泥沙起动流速公式存在的不足并对此作了修正[8]。黄永健等（1997）结合灌渠的结构特点及水流泥沙特性, 对一维非耦合、非均匀、不平衡输沙基本方程组进行了合理的简化, 根据实体模型试验及实测资料, 确定了自流与提灌时干、支渠分流的不同分沙分水比并应用于黄河下游的输沙计算[9]。迄今为止的输沙计算方法普遍存在两个不足，一是经验型公式（模型）较多，其应用范围常被限制在研究开展的流域或河段；二是计算精度普遍较低，对输沙过程的模拟精度不高。本研究探讨将MacCormack方法用于推移质的输沙计算，为河流输沙，河床演变以及河道疏浚等提供科学的计算方法。

2 MacCormack方法

MacCormack方法由Robert W. MacCormack在1969年提出，是计算流体力学中用来求解双曲形偏微分方程数值解的一个普遍离散方法，为二阶有限差分方法。该计算方法非常简洁，易于理解和实现编成计算，应用MacCormack方法计算时包括两个阶段：预测阶段和修正阶段，该方法非常适合线性方程的求解。由于采用二阶差分对预测阶段和修正阶段进行计算，数值计算结果的精度较高。对于非线性方程，应用MacCormack方法也可得到较好的结果。

3 算法建立

3.1 基础方程式

用于径流输沙过程计算的基础方程式包括地表径流连续方程式，运动方程式以及输沙量计算公式[10]。

·地表径流连续方程式（一维稳定流）

（1）

·运动方程式

（2）

·输沙量计算公式（Ashida / Michiwue 经验公式）

（3）

式中，Δt为计算的时间步长，1s；Δx为河流纵向计算的单位步长，50 m；b为河道宽度，m；Q为流量；g为重力加速度，m/s2；h为水深，m；ib为河道平均坡度；ie为水力坡度；qbi为单位时间内单宽河道输沙量，m3/s；s为沙在水中的比重，kg/m3；d为流沙的平均粒径，0.5 mm，为推移质；τ*c为界限推移力，kg；τ*为推移力，kg；u*c为界限摩擦速度，m/s；u*为摩擦速度，m/s。

3.2 有限差分

计算时，需要依据地表径流连续方程式及运动方程式耦合输沙计算在时间上离散化，即有限差分。因为计算依赖于一定的初始条件和边界条件，所以预测阶段采用后退差分法，修正阶段采用中间差分法[11,12]。

·预测阶段差分式

（4）

·修正阶段差分式

（5）

式中，i为计算的栅格编号，n为计算的时间编号，A为过流面积，m2; kv为人工粘性系数，0.00001；其它因子同上。

3.3 计算程序编译

实现数值计算的程序采用计算机高级语言Fortran编译，此部分内容在此略去。

4 计算事例

4.1 计算条件

图1 计算河道长度及形状

计算河道的形式如图1所示，河口前长度为2 km，河床坡度ib为0.001。自河口以下计算长度为1k m，水流自河口开始以与河道中心线成20o角向两侧扩散，河口以下坡度为0.01， 河口的标高为0。河道年均流量为 50 m3/s，计算开始5年后发生了100年一遇的洪水，其设计流量为2500 m3/s，洪峰持续时间为5小时。求解问题为5年后发生洪水前河床的形态，洪水发生后每隔1小时河床的演变情况。

4.2 结果及讨论

利用上述的数值计算方法和题设条件进行计算，根据求解问题的要求，当年平均流量为50m3/s时 5年后的河床形状、水面线以及初始状态的河床形式如图2所示。

图2 河床的初始形态与5年后形态的对比

洪水发生后，每隔1小时河床的形态与水面线的数值计算结果如图3所示。

图3洪水发生后每1小时间隔的河床及水面线形状

由图2可知，因输沙的影响，河床的初始形态与5年后的形态明显不同，特别是河口部发生了明显的堆沙现象。

图3描述了洪水发生后，水面线及河床变化的过程，河床形态以及水面线在河口部以前呈现相对稳定的状态，没有发生明显的变化，即河床演变在河口部上游区域表现并不明显。在河口部，洪水发生后导致推移质随时间逐渐前移，并在河口部堆积，其泥沙的堆积区域主要在距河口700 m范围之内，而后泥沙堆积的程度逐渐变小，在距河口1000m之处几乎没有发生泥沙堆积现象。

利用数值计算计算的结果，可以推求输入河口的泥沙量，即河口的堆沙量。其计算方法是比较河口以下发生洪水后的河床高程与计算初始时河床各段的高程（Δx=50 m），由二者的差值结合河道的宽度可近似计算出输沙量。从而为河口部泥沙疏浚提供基础数据。

5 结论

研究实现了MacCormack方法与推移质输沙计算耦合，为河流推移质计算提供了一种实用的方法。

计算结果可以很好地模拟河道推移质移动的过程，利用MacCormack方法进行河流输沙计算可以得到准确度较高的结果。

研究可为河流输沙量的准确推求，河床演变以及河道疏浚等提供数值计算方法。