高能河流二维输沙模拟在渭河上的应用

董礼玮

（西安水文水资源勘测中心,陕西 西安 710000 ）

1 引言

地下水是农村和城市地区最重要的水资源之一。由于淡水资源与河流在数量和质量上的强烈交换,淡水资源极其脆弱。随着与气候变化相关的水资源的减少和需求的增长,大多数国家将地下水和地表水视为单一资源[1-2]。事实上,这两个实体是相互关联的,根据其数量或质量,它们之间关系密切。河流形态动力学是对控制和维持河流环境的物理过程的重要研究。因此,了解这些复杂的形态变化是了解河流和含水层之间交换区的关键问题。河流汇流对于河流系统中的水和沉积物的路由以及形态演变非常重要[3-4]。作为河网连接的河道中的分流现象通过这些交叉点相互影响,比单个河道中的分流现象复杂得多。对于沙质河流而言,在高含沙汇合洪水期间,其形态动力学在时间和空间域上都会发生显著变化,这意味着水流运动和泥沙输移之间存在着积极的相互作用,这给河流模型的开发带来了更大的挑战。因此,利用水动力模型模拟高浓度汇流洪水仍然是一项困难的任务。

2 研究区概况

2.1 渭河下游

渭河下游是黄河最大的支流。研究段长134.4 km,宽6 km~10 km,河道纵坡0.3‰~0.6‰,北洛河在潼关上游12 km~20.5 km处注入渭河。三门峡水库蓄水后10 年间,渭河下游河道曲流频繁,北罗河入口位置发生突变。冲积含水层是当地居民的主要饮用水资源之一。事实上,最近几年的统计结果提到,有5000 万m3的水是从这种非承压含水层中抽取的[5]。几十年以来,人类活动（沉积物提取、城市化等。）影响着地下水位的变化。

2.2 渭河下游形态动力学的关键问题

渭河下游的农业发展、城市化和沉积物提取导致了河流形态变化。事实上,为改善渭河下游农业发展等,修建了许多水利工程。修筑堤坝是为了保护人口和农业用地免受洪水侵袭；为降低河床,从河床中提取了1500万m3（相当于150年的自然沉积物供应）沉积物。此外,还建造了拦水堰,通过地表水下渗来增加地下水位。河床从150 m延伸到280 m,这导致了水流速度的增加。因此,在下游沿河许多地方可以观察到河岸侵蚀不断发生,这些物理过程影响了河流泥沙和水流之间的运移情况。

除了人为活动引起的河道形态变化外,河床的物理动力受洪水期影响,导致沉积物重组。为了解形态动态及其对河流含水层交换区的影响,必须建立数值模型。本文介绍了实现数值模型的方法和用于模拟渭河下游河道形态变化。

3 数据和方法

河流演变的建模需要计算模块来模拟流体动力学和床层演化。如果以适当的方式管理开发过程,2D模型可以将演变过程进行可视化。

3.1 数据

2009年~2013年期间渭河下游河势变化频繁,不仅河流形态发生了变化,洪水频发,导致渭河流域径流量呈枯丰交替。利用2011年和2013年的两个数字高程模型对渭河下游的形态动态进行评价,这种高分辨率和全维度的河道地形模型是开发表征河流特征和理解水流挟沙运动机制的新方法,为河流输沙数值模拟提供了强有力的信息。2011年洪水事件引起的形态变化可以通过比较两个地面模型来识别（图1）。虽然这种方法在几何角度方面有很多优点,但是受数据的准确性影响,这种方法也有一定的局限性。

图1 2011年洪水事件中识别沉降和侵蚀的比较

3.2 流体动力学模型

MIKE 21 FM的HD模型,通过求解质量和动量的二维浅水方程,进行了流体动力学模拟。为了使用有限体积数值方法求解这些方程,使用非结构化网格来定义研究区域的地形。模拟在单层流体中产生不稳定的二维流动（垂直均匀）。流量和水位变化由Saint-Venant方程（垂直方向上积分的质量和动量守恒）描述,如下所示：

式中：h是水深,m；w=R-I是净流量,m/s；z是表面高度,m；（p,q）=（hu,hv）分别是x和y方向的通量密度,m2/s；g是重力加速度,m/s2；pa是大气压,kg/（m·s2）；pw是水密度,kg/m3；τxx,τyy,τxy是有效剪切应力的组成部分。

3.3 泥沙输送模型

使用MIKE 21 FM中的输沙模块（ST）进行了床层演变的模拟,该模块计算沉积物输送能力,河床层水平输沙的初始速率以及由于水流运动引起的非粘性沉积物的形态变化。沉积物输送计算基于流体动力学条件和沉积物性质。ST模块考虑了河床层负荷和悬浮负荷。河床层负荷由剪切应力或每单位流量控制,并与流动瞬间反应。悬浮负荷是与流量相比,运输中的相位滞后。河床荷载不需要对流色散模拟,但必须考虑两个重要因素：河床剪切应力方向的偏差和倾斜河床的影响[5]；另一方面,流体中非粘性悬浮沉积物的建模可以通过体积沉积物浓度的运输方程来描述。

式中：Sbl为推移质,m3/s；Ssl是悬移质,m3/s；T是取决于临界摩擦速度和有效摩擦速度的无量纲传输阶段；D*是无量纲粒子参数；S是沉积物的相对密度；Ca体积床浓度,m；f基于床层荷载的悬浮荷载修正系数；V是速度,m/s。

确定水平变化的关键参数是元素单元中心的床层水平变化率。河床的这种演化是用沉积物连续性方程Exner方程得到的：

式中：n是河床的孔隙度；z是河床的高度,m3；Sx是河床的推移质或者总的推移质沿着x方向,m2/s ；Sy是河床的推移质或者总的推移质沿着y方向,m2/s；ΔS是泥沙的汇或者源。

根据估计的床面变化率,通过形态学模拟对床面进行连续更新。泥沙输送部分的目的是比较建立一个模型的合适方法,以MIKE 21 FM和ST模块能够代表渭河下游河道的形态变化,以及用于求解河床形态变化的方程。

4 结果和讨论

4.1 流体动力学结果

首先,比较不同几何形状对流体动力学模拟的影响,根据网格分辨率和作业管理的需要,讨论了5 种分辨率的网格模拟。为了得到一个高效的MIKE 21 FM 模型,对同一个洪水事件进行了多个网格的模拟,研究区域采用不同分辨率（5 m、10 m、15 m、20 m、25 m）的三角形离散化表示。然后将三角形网格与四边形网格进行了比较。最后,利用 villemonte 堰流公式实现了水工建筑物的施工计算,这是一个经验公式[6]。本文选择2015年10月3日~2015年10月6日一系列的洪水事件进行模型的测试。模型与观测资料的比较（图2）表明10 m的分辨率是最准确的。采用三角形离散化方法模拟溢流堰的流动更为有效（图3）。模型中溢流堰结构化的实现需要许多参数,并不能真正代表堰效应（图4）。最后,将三角形网格,结构的三角形网格和四边形网格的模拟速度进行比较,得出三角形网格以及四边形网格模拟效果要优于结构三角形网格模拟。

图2 不同分辨率的模型与观测资料的水深比较

图3 模型中采用三角形网格与四边形网格计算比较

图4 模型中采用三角形网格与结构三角形网格计算比较

图5 各网格的速度比较

4.2 河床形态演变

渭河存在的主要问题之一是下游淤积。泥沙输送模型考虑了推移质和悬移质的影响。因此,2011年洪水事件的模拟再现了包括堰在内的河床水位。流体动力学研究结果表明,结构的实施不同程度地影响流场。针对堰体对地貌的影响,通过对三种模型的比较,找到解决方案: 一是模型以均匀的颗粒尺寸定义,二是模型以结构化网格划分实现,三是模型以非均匀的颗粒尺寸,可定义为阻止破坏堰体的推移质。三种模型的河床变化比较（图6）表明,对于渭河河床模拟的这一部分来说,结构的实施是模拟河床物理过程的最合适的方法。河床演变的计算是基于水动力学和泥沙特性。在这里,只有沉积物的特性不足以描述河床形态变化。

图6 各模型模拟的床层变化比较

5 结语

渭河由于频繁的高含沙洪水,河道演变十分活跃。基于水动力和形态过程耦合在高浓度合流洪水模拟中的重要性,本文提出2D水力模型对河流内泥沙输移过程进行模拟。利用水动力学模型的离散化方法,建立了泥沙输送模型。模拟结果表明,三角形网格（10 m分辨率）虽然是描述堰上流场的最佳方法,但不适合模拟堰下的河床演变。因此,它需要在实施模拟的结构的基础上提高模拟性能,同时结合相关河流输沙经验公式法,应该着重分析该模型使用不同泥沙输送公式的灵敏度,并将这种方法应用河于道的蜿蜒部分。