清水沟水库溃坝对主河道行洪过程影响数值模拟研究

马利平， 侯精明， 刘昌军， 李桂伊， 高 波， 蔚桂玲

(1.西安理工大学 省部共建西北旱区生态水利国家重点实验室， 陕西 西安 710048； 2.中国水利水电科学研究院， 北京 100038； 3.西安地质调查中心， 陕西 西安 710054; 4.青海省海北藏族自治州水利勘测设计室， 青海 海北 812200)

1 研究背景

小型水库是我国解决农村居民生活用水、防止水土流失、发展农村经济状况的重要工程措施。20世纪50-70年代，我国建成小型水库75 646座，其中近95%为土石坝[1]。这些土石坝普遍存在工程质量差、防洪标准低、建成年代久远、水库淤积严重等问题。一旦失事溃决，将给下游带来巨大的灾难[2]。因此高精度的溃坝洪水模拟预警就显得尤为重要。

在溃坝洪水演进模拟研究中，模型对于复杂地形下干湿边界的处理[3-4]、界面通量的计算能力直接决定着模型的稳定性和精度，因此，处理好干湿边界及界面通量的计算是溃坝模拟研究的重点问题[5]。针对这一问题，吴钢峰等[6]基于结构网格，采用有限体积法构建了二维水动力学模型，采用中心迎风格式求解界面通量，模拟了溃坝洪水在复杂实际地形条件下的流动过程。夏军强等[7]在无结构网格下建立了采用有限体积法的二维水动力模型，利用Roe格式的Riemann解计算界面水流通量，研究了不同最小水深对洪水演进的影响；向波等[8]基于三角网格建立了求解浅水波方程的有限体积模型，很好地捕捉了溃坝波的前进。WANG等[9]提出一种多线截面图法，研究了在不规则断面上溃坝波的解析解；然而，在现有的众多溃坝洪水模型中，并未将负水深带来的稳定性问题和求解界面通量的精度问题达到较好的和谐统一。HOU等[10-11]采用有限体积法建立了求解二维浅水方程的水动力模型，准确模拟了Malpasset溃坝事故，但尚未应用其模拟支沟对主河道的影响。在众多溃坝事故中，有不少水库位于支沟，支沟溃坝洪水对主河道洪水肯定有影响，但该影响如何的研究却很少。

为研究支沟溃坝对主河道洪水过程的影响，本文以发生于2017年7月陕西省榆林市清水沟水库溃坝为研究对象，采用了集成HLLC近似黎曼求解器的二维水动力模型，模拟了位于支沟的清水沟水库溃坝对主河道大理河洪水过程的影响，并系统分析了主河道上下游行洪过程对支沟溃坝的响应规律。

2 模型简介

2.1 控制方程

本文应用二维非线性浅水方程来模拟溃坝洪水的演进过程，其矢量形式[12]如下：

(1)

(2)

式中：q为变量矢量，包括水深h，两个方向的单宽流量qx和qy；g为重力加速度;u和v分别为x、y方向的流速；F和G分别为x、y方向的通量矢量；S为源项矢量；zb为河床底面高程;Cf=gn2/h1/3为谢才系数;n为曼宁系数。

2.2 数值方法

本模型采用Godunov格式的有限体积法求解二维浅水方程。单元界面通量应用HLLC格式的近似黎曼求解器求解[13]。通过静水重构来修正干湿边界处负水深。底坡源项使用模型作者提出的底坡通量法处理。摩阻源项的计算使用半隐式法来提高稳定性[14]。采用二阶显式Runge Kutta方法进行时间步长的推进[15]。从而构造具有二阶时空精度的MUSCL型格式，有效解决复杂地形干湿界面处的负水深和伪高流速等现实中不存在的物理现象所造成的计算失稳和物质动量的不守恒。因溃坝洪水演进过程模拟一般尺度较大，为提高计算效率，模型采用GPU并行计算技术实现高速运算[16]。

2.3 黎曼求解器

Godunov型Riemann求解器普遍用于计算非线性双曲线方程(如SWEs)的通量。计算公式如下[15]：

(3)

(4)

式中：vL和vR为左右切速度分量。SL、SM和SR分别为左、中、右波速，其计算公式如下：

(5)

(6)

(7)

0.25(UL-UR)]2

(8)

利用黎曼近似求解器求得F*：

(9)

3 模型构建

3.1 研究区域概况

清水沟水库位于陕西省榆林市子州县清水沟与大理河汇合口上游约100 m处(图1)，是子洲县供水水源工程，工程以上流域面积5.97 km2。子洲县城所处区域地下水均为苦咸水，水质较差，不适于居民生产生活用水。该工程为子洲县城主要供水水源，工程运行方式为在丰水期将大理河的河水抽蓄至水库存蓄，再由放水管和输水管道向县城供水。该工程主要包括大坝、抽水工程、放水工程。

3.2 模型设置

模型的输入资料为地形数据、入流数据等。根据溃坝发生后的实际地形，对于支流溃坝段采取40 m宽的入流边界，对于干流大理河段采取70 m宽的入流边界，下游边界为自由出流的开边界，其余边界定义为闭边界。

图1 研究区域水系及位置图

3.2.1 地形资料 溃坝通常会伴随河道冲刷并引起地形改变，但由于这一过程极其复杂，很难模拟，目前国内外对此问题的研究较少，大多数学者都采用溃坝前的地形或溃坝后的地形进行溃坝洪水的模拟计算[7,10,17]。考虑到本文所研究的是清水沟支沟溃坝对大理河主沟行洪过程的影响，冲刷引起的地形改变对溃坝大流量的影响不大，故采用溃坝后的地形进行模拟研究，采用溃坝前的地形进行模型验证。

研究区域面积为0.274 km2，为此次溃坝所能影响的主要区域。溃坝发生后，中国水利水电科学院研究团队采用三维激光扫描仪，对溃坝区域进行了三维激光扫描，具体范围从大理河口到上游水库回水区的末端，获取了该区域内高分辨率的全彩色LiDAR数据，经处理得到精度为1 m网格数为274 320的数字地形高程数据(图2)，图中溃坝处至大理河交汇处长135 m，上游河道长373 m，下游河道长294 m。

3.2.2 入流资料 2017年7月26日位于大理河子州县城上游的清水沟水库发生溃坝，库水泄入大理河。根据陕西省防汛抗旱总指挥部办公室报告，在溃坝发生时大理河子洲段流量约为800 m3/s，故给予干流大理河段800 m3/s的恒定流量作为入流条件，持续75 min。根据陕西省防汛简报、大坝现场工作人员拍照数据以及陈祖煜等研发的溃坝洪水分析软件(DB-IWHR2014)[18]综合得到了溃坝洪水过程线，将此溃坝洪水过程线给予支流溃坝段作为入流条件，水库从溃坝开始至水库库容泄空共持续35 min。溃坝洪水过程线如图3所示。

图2 研究区域数字地形高程及代表断面位置图

3.2.3 参数设置 本次模拟共模拟75 min，在未发生溃坝之前，给予干流大理河段800 m3/s的恒定流量，持续30 min。待大理河段流量稳定，给予支流溃坝段溃坝流量，持续35 min。直到水库库容泄空，再模拟10 min等待研究河段流量稳定。根据现场实际情况选取曼宁系数为0.02[19]，克朗数为0.5。

3.3 模型验证

文献[10,20-21]中Hou Jingming等应用此模型模拟研究了Malpasset溃坝洪水演进过程，模拟结果与下游各实测水位吻合良好，验证了此模型的可靠性。本次模拟根据中国水利水电科学研究院后期现场调研结果显示，洪水最大淹没高度为955.3 m，即图4(a)中洪痕位置。应用本模型模拟时，由于溃坝引起坝后地形的冲刷，故采用溃坝前的地形进行模型验证，选取支沟下游与图4(a)中洪痕附近位置的水位(图4(b))，得到此处的水位-时间曲线，如图4(c)所示，最大水位为955m。故模拟结果与实际洪痕较为一致。

4 模拟结果分析

4.1 溃坝洪水过程模拟

本次溃坝洪水过程共模拟75 min，在前30 min给予大理河段800 m3/s的恒定流量，待河道流量稳定(图5(a))。在30 min后，溃坝洪水开始从支沟向大理河演进，经过100 s，溃坝洪水到达与大理河交汇口处(图5(b))，第38 min溃坝洪峰出现(图5(c))，之后溃坝流量逐渐减小至溃坝结束，总历时35 min。图6为洪水演进过程的速度矢量图，图中箭头大小代表该处流速的大小，箭头方向为该处流速的方向。图7为溃口处的溃坝波演进过程。

图4 下游洪痕位置图及洪痕断面水位图

图5 溃坝洪水过程模拟

图6 溃坝洪水速度矢量图

图7 溃坝波演进过程模拟

4.2 上游断面洪水过程

选取大理河上游1号、2号、3号断面作为代表研究断面(图2)。

对各断面的水位，由于溃坝洪水到达大理河干流时有一部分洪水向上游演进，故水位曲线有先逐渐升至最高再缓慢下降的趋势，模拟结果发现，上游各断面水位抬升的范围在0.392～0.404 m之间。对于上游各断面的流量，由于有一部分溃坝洪水向上游演进，故流量曲线有先下降后上升再下降至稳定的趋势，模拟结果发现，此种趋势随着距溃坝洪水汇流处越远越不明显，即对于干流上游来说，距溃坝处越远，流量对于上游的影响越小。

4.3 下游断面洪水过程

选取大理河下游4号、5号断面作为代表研究断面(图2)，对于各代表断面的水位，由于溃坝洪水到达干流大理河时大部分洪水向下游演进，故水位曲线统一有先逐渐升至最高再缓慢下降的趋势，模拟发现，下游各断面水位抬升的范围在0.325～0.55 m之间。对于下游各断面的流量，由于大部分溃坝洪水向下游演进，故流量曲线有先迅速上升再缓慢下降至稳定的趋势，模拟发现，下游各断面处最大洪峰流量为982 m3/s，相比简单的叠加，洪峰流量减少了18 m3/s。

由4.1节可知该二维水动力模型能较为准确地模拟溃坝洪水的演进过程。其中溃坝洪水速度矢量图及溃坝波演进过程图可较为直观地看出溃坝洪水泄流而下汇入大理河的情况。上游的水深流量图表明了支沟溃坝在主河道上游产生的壅水现象及流量减小的事实，同时结合下游的水深流量过程，可以发现支沟溃坝对主河道上下游都有影响，其中最为关心的下游洪峰流量并不是简单的叠加，而与上下游关系密切。

图8 大理河上游断面水深、流量过程线

图9 大理河下游断面水深、流量过程线

5 结 论

研究支沟溃坝洪水在主河道的行洪过程对于干流上下游的防灾减灾具有非常重要的意义，本文应用集成HLLC近似黎曼求解器的二维水动力模型，模拟了清水沟水库溃坝对主河道大理河洪水过程的影响，得出以下结论：

(1)本模型可用来模拟溃坝洪水的演进过程，模拟研究发现，此次子州清水沟水库溃坝并未对下游造成明显的追加破坏，这与实测资料是一致的。原因是溃坝发生时，溃坝洪峰流量200 m3/s直接汇入大理河，但此时大理河的流量只有800 m3/s，先前3 150 m3/s的洪峰流量已经通过。

(2)溃坝洪水对大理河上下游均有影响，在上游处会产生明显的壅水现象，导致水面增加了0.392～0.404 m，同时流量具有先下降后上升再下降至稳定的趋势；其中最为关心的下游洪峰具有先增加再下降至稳定的过程，洪峰的传播速度为2.06 m/s，在断面5处的洪峰流量为982 m3/s，相比简单的叠加，洪峰流量减少了18 m3/s，原因是溃坝波向上游演进并最后由于重力作用下移，最大洪峰不是简单的叠加，而与上下游关系密切。该成果对于防治支沟溃坝在主河道形成洪水具有一定的指导意义。