不同频率洪水遭遇下溃坝洪水演进的二维模拟

冯民权，杨 志

(西安理工大学 陕西省西北旱区生态水利工程重点实验室，陕西 西安710048)

溃坝是一种低概率、高危害的社会致灾因素，大坝一旦失事其危害性极大。为分析溃坝洪水风险、评估灾害损失以及制定地区应急预案，研究溃坝洪水演进十分必要。溃坝问题研究方法始于物理模型试验[1-2]，但该方法对于存在大量边界形状复杂区段的实际工程，往往只能给出某一区段总流的参数；而通过一维[3]、二维[4]数值模拟则能给出相关流场的具体信息，因此数值模拟方法渐已成为溃坝流体力学研究的强有力手段。如Baghlani[5]基于笛卡尔网格有限体积法，研究了一维及二维非对称缺口、环状溃坝、洪泛区设置障碍等溃坝模拟问题，均得到了满意的结果。在国内，周孝德等[6]在黑河金盆水库洪水淹没范围及撤离方案研究报告中，对不同频率工况下黑河洪水的演进过程进行了一维、二维模拟计算与分析。目前，尽管溃坝数值模型研究已取得了长足的进展，但仍存在有待深入探讨的问题。第一，采用结构网格法[7]模拟溃坝较多，而无网格法如粒子法[8]、SPH法[9]等应用较少。然而，溃坝发生处的地形往往是复杂的不规则图形，而结构网格法对地形限制较大，无网格法尚不成熟，采用非结构网格更加符合地形需要，也具有更高的精度。第二，就黑河金盆水库而言，该水库位于秦岭北麓山区渭河支流上，水库下游淹没影响区一方面具有大面积的平原地形特点，需采用平面二维非恒定流数学模型模拟溃坝洪水的演进过程；另一方面，计算区域涉及部分山区地形，其高程跨度大(384.7～975.5 m)，溃坝水流水位变幅大，需及时调整边界。第三，不同频率入库洪水相遭遇的调洪演算及其漫顶导致溃坝洪水的演进模拟有待进一步研究[10]。针对以上问题，本研究参考前人研究[6,11]成果，试图建立合适大小及密度的非结构网格模拟金盆水库溃坝洪水的演进，并根据实际地形调整边界属性，确保模型计算过程的稳定；优化各部分计算方法，确定溃坝模拟参数，采用MIKE软件进行计算，建立黑河金盆水库下游10 000年一遇入库洪水漫顶失事导致的溃坝洪水在下游演进的二维模型，并在此基础上分析100年一遇山区暴雨洪水延时遭遇10 000年一遇入库洪水的淹没水深、范围及流速变化，旨在为黑河金盆水库的安全运行和调洪决策提供参考，并为洪水风险分析和灾害损失评估提供依据。

1 数学模型

1.1 控制方程

由于溃坝洪水运动的水平尺度远大于垂直尺度，水深、流速等水力参数沿垂直方向的变化较水平方向的变化小得多，因而可将三维流动控制方程沿水深积分，并取平均水深，从而得到二维浅水水流运动方程。

水流连续方程：

(1)

式中：x为水平方向的空间坐标，m；y为垂直方向的空间坐标，m；Z为点(x,y)处的水位，m；t为时间，s；h为点(x,y)处的水深，m；u为点(x,y)水平方向的流速分量，m/s；v为点(x,y)垂直方向的流速分量，m/s。

水流运动方程：

(2)

(3)

1.2 方程的离散

针对溃坝洪水数值计算中面临的不规则边界和复杂地形等问题，建立非结构三角形网格下的二维水流模型，将任意三角形作为计算单元，采用单元中心的有限体积法求解方程，单元网格的分布如图1所示。其中，点M代表各项结果变量位于单元中心，F为跨边界通量，分别垂直于所在三角形单元中相应的各条边。

图1 非结构三角形网格及其单元的分布

对浅水水流运动方程进行归一化，则有：

(4)

对方程(4)在每一个单元上进行积分，经进一步简化后，得：

(5)

式中：Ui为控制单元i的守恒变量；Ai为控制单元i的面积，m2；N为控制单元的边的个数；F为通过第j条边的通量；n为通过第j条边的法向量；ΔΓj为控制单元各边长度；Si为控制单元i的源项。

用Roe’s近似黎曼解法求解跨边界通量F，用二阶TVD格式的限制器函数避免数值振荡。

1.3 边界条件

根据黑河金盆水库实际情况及已有资料，二维模拟入流边界条件采用坝址断面的流量变化过程Q=Q(t)，其y方向流速假设为零。出流边界采用开边界条件。

针对农村初中古诗词教育轻视诵读教学、照搬参考书、忽视激发学生兴趣、教学手段单一等问题，依据“教师为主导、学生为主体”的理念，通过课内进行情境教学改革，课外进行初中语文古诗词读书分享学习模式探索以及“网络多媒体+古诗词”的学习方式研究，形成“课内+课外+网络多媒体”三位一体相结合的教学模式。提高了学生对诗人情感的理解能力，培养了学生的古诗词鉴赏能力，提高了学生的人文素养和审美情趣，引导学生认识中华文化的丰厚博大，吸收民族文化智慧，激发学生的想象力和创作潜能，促进个性发展，丰富精神世界。

2 计算条件

2.1 计算区域

黑河金盆水库位于西安市周至县城南黑河干流峪口以上1.5 km处，东距西安市约86 km，北距周至县城约14 km，坝址处控制黑河流域面积的65.6%。水库设计正常高水位为594.0 m，总库容2.0亿m3、有效库容1.77亿m3，溃坝计算类型按照洪水漫顶失事拟定，即当静水位等于600 m时发生溃决。计算区域涉及黑河金盆水库坝址至黑河入渭口及渭河部分流域，总面积为1 066 km2，沿程地形及A、B、C、D、E 5个观测点位置如图2所示。区域内采用三角形网格划分，并通过网格密度控制优化处理突变边界，最终确定16 199个网格，每个网格面积约为0.065 8 km2(图3)。

2.2 糙率的选取

由于黑河金盆水库地形复杂，糙率值难以准确估计，而糙率选取的恰当与否对计算成果影响较大。糙率的取值范围一般为0.025～0.045[12]，参照文献[11]，最终拟定二维计算的河床糙率值为0.031。

图2 黑河金盆水库计算区域地形图

图3 黑河金盆水库计算区域网格剖分图

2.3 工况设计

设定10 000年一遇入库洪水漫顶导致的溃坝洪水为工况1。据现有水文资料记载，1898年秦岭山区大面积暴雨形成渭河中、下游洪水，同期黑河支流调查流量达3 620 m3/s，为该河段近百年来最大一次洪水，在黑河历史洪水中比较具有代表性，故选100年一遇入库洪水延时遭遇10 000年一遇入库洪水为工况2。在上述2种工况下黑河水库入库洪水的流量过程线如图4所示。

图4 工况1和工况2下入库洪水的流量过程线

2.4 水库调度方式

(1)当水库水位小于汛限水位591 m时，下泄流量为0。

(2)当水库水位在汛限水位591 m至防洪高水位594 m之间时，如果来水流量小于下游允许泄流量1 800 m3/s，按总下泄流量等于来水流量调洪；如果来水流量超过下游允许泄流量1 800 m3/s，控制总下泄流量等于1 800 m3/s。

(3)当水库水位大于防洪高水位594 m时，闸门全部开启。

2.5 边界条件

工况1：由金盆水库10 000年一遇入库洪水的流量过程线(图4)可知，此时最大入库流量为7 800 m3/s。经过调洪演算得到水库下泄流量过程线见图5。从图5可以看出，从t=0时刻开始调洪，经过16.4 h，水库水位达600 m，大坝溃决。之后峰前流量经过1.67 h上升阶段，于t=18.1 h至决口口门达到最大(表1)。峰后流量过程线概化为4次抛物线型进行计算，结果见图6。

图5 工况1洪水调洪演算的下泄流量过程线

表1 工况1峰前流量过程线

综合上述计算结果，整理得到坝址处下泄流量见图7，以此作为二维计算时入流的边界条件。 工况2：如图4所示工况2的入库洪水流量过程线，从t=0时刻开始，100年一遇洪水入库，经调洪演算后下泄，该过程持续至t=23 h，此时10 000年一遇洪水入库共同参与调洪下泄，该过程持续至t=37.5 h，即刻水位升至600 m，大坝溃决。采用工况1所述方法进行计算，综合整理得到入库流量叠加过程及叠加后的下泄流量过程见图7。

全部计算区域出流边界设定为开边界。在黑河入渭口加入限定条件，计算所得流量水位关系(Q=Q(z))见图8。

图7 工况1、工况2下泄流量过程线

3 计算结果与分析

取黑河金盆水库坝址处为上游进口断面，以黑河入渭口为下游出口断面，其余全部设置为开边界，洪水标准为10 000年一遇，经过水库调洪演算，当库水位与坝顶齐平时认为发生溃坝事件，即溃坝模拟计算开始。选取计算步长Δt为8 s，计算总时长为73 h，对10 000年一遇溃坝洪水的演进过程进行二维模拟分析。

3.1 工况1计算结果

图9为溃坝洪水演进1，2和3 h时的水深等值线图。

图9 工况1溃坝后1，2，3 h时洪水淹没水深的分布

从图9可以看出，溃坝发生一段时间后，洪水最大淹没水深很快出现下降趋势，而同水深值的淹没范围显著增加，说明溃坝洪水不仅具有急骤破坏性，同时具有大范围破坏性，因此对溃坝洪水进行实时演进模拟意义重大。图10为溃坝洪水演进1和2 h时的流场图，可见淹没范围内的水流流向呈放射状扩展规律，沿水库向下游河道演进，流速呈中心较急，两边趋缓分布。

图10 工况1条件下溃坝发生1，2 h后洪水流速的分布

图11直观地表现了溃坝洪水在大坝下游黑河淹没区A、B、C、D、E 5个观测点的演进过程。由图11各观测点水深变化过程和流速变化过程可知，溃坝洪水到达A、B、C、D、E的时间分别为模拟初始时刻后的16.45，16.68，16.97，17.58和18.01 h，最大淹没水深分别为6.18，3.31，3.10，4.39和3.92 m，最大流速分别为8.68，4.47，4.46，2.28和2.35 m/s。可见，随着下游地区与坝址处距离的增大，洪水到达时间变长，洪水最大流速变小，溃坝洪水急骤性损害能力逐渐降低。

图11 工况1条件下下游各观测点水深和流速的变化过程

3.2 工况2计算结果及与工况1的对比分析

图12为计算区域在工况2条件下溃坝不同时刻的水深分布。从图12可以看出，在同一时刻，工况2比工况1淹没范围更广，淹没水深更大，差别主要表现在下游区域；不同时刻，溃坝洪水演进规律与工况1结果分析中所述相符。

图13为工况2条件下黑河金盆水库下游A、B、C、D、E 5个观测点水深和流速的变化过程。图13显示，溃坝洪水到达A、B、C、D、E的时间分别为模拟初始时刻后的16.02，16.18，16.51，16.82和 17.17 h，最大淹没水深分别为6.91，3.66，4.01，4.41 和 4.19 m，最大流速分别为9.40，4.81，4.46，2.36和2.46 m/s。与工况1相比，工况2淹没水深变化较大，且会使溃坝最大流量时刻提前出现；流速分布的变化较小，但同一时刻、同一观测点的流速均大于工况1。

图12 工况2条件下溃坝1，2，3 h后水深的分布

图13 工况2条件下下游各观测点水深和流速的变化过程

4 结 论

1)本研究针对黑河下游地区高程跨度很大的地形特点，建立非结构网格系统，调整三角网格的大小和密度以及边界属性，确定溃坝计算参数。对控制方程进行离散及简化，通过Roe’s近似黎曼解法对之求解，从而建立了黑河金盆水库下游溃坝洪水演进的二维模型。

2)绘制了模拟工况1条件下黑河水库漫顶溃坝洪水不同时刻的水深分布图及流速矢量图，得到溃坝洪水随时间的变化过程，可知随溃坝发生时间的延长，洪水的急骤破坏性有所降低，但破坏范围激增，水流流速表现为中心急、两边缓的放射状扩展规律。在拟定的工况2条件下，与工况1相比其淹没范围更广，淹没水深显著增加，溃坝洪水最大流量出现时刻提前；流速分布变化不大，但工况2下同一时刻、同一观测点的流速均大于工况1。

3) 利用本研究的计算成果可以确定2种模拟工况下黑河下游地区不同观测点的洪水到达时间和最大淹没水深，为下游淹没区域的防洪规划、预警决策等提供了科学依据。

[参考文献]

[1] Kocaman S,Cagatay H O.The effect of lateral channel contraction on dam break flows:Laboratory experiment [J].Journal of Hydrology,2012,432/433:145-153.

[2] Carrivick J L,Jones R,Keevil G.Experimental insights on geomorphological processes within dam break outburst floods [J].Journal of Hydrology,2011,408:153-163.

[3] Chang T J,Kao H M,Chang K H,et al.Numerical simulation of shallow-water dam break flows in open channels using smoothed particle hydrodynamics [J].Journal of Hydrology,2011,408:78-90.

[4] Singh J,Altinakar M S,Ding Y.Two-dimensional numerical modeling of dam-break flows over natural terrain using a central explicit scheme [J].Advances in Water Resources,2011,34:1366-1375.

[5] Baghlani A.Simulation of dam-break problem by a robust flux-vector splitting approach in Cartesian grid [J].Scientia Iranica,2011,18(5):1061-1068.

[6] 周孝德,冯民权,闵 涛,等.黑河金盆水库洪水淹没范围及撤离方案研究报告 [R].西安:西安理工大学,2006.

Zhou X D,Feng M Q,Min T,et al.The report of the submerged area and evacuation plan of Heihe Jinpen reservoir flood [R].Xi’an:Xi’an University of Technology,2006.(in Chinese)

[7] 沈 洋,王佳妮.基于MIKE软件的溃坝洪水数值模拟 [J].水电能源科学,2012,30(6):56-58.

Shen Y,Wang J N.Numerical simulation of dam break flood based on MIKE software [J].Water Resources and Power,2012,30(6):56-58.(in Chinese)

[8] Shakibaeinia A,Jin Y C.A mesh-free particle model for simulation of mobile-bed dam break [J].Advances in Water Resources,2011,34(6):794-807.

[9] 缪吉伦,赵万星,黄成林.SPH 法模拟立面二维溃坝流动应用研究 [J].重庆交通大学学报：自然科学版,2012,31(1):121-123.

Miao J L,Zhao W X,Huang C L.Application of vertical Two-Dimensional dam break simulated by SPH method [J].Journal of Chongqing Jiaotong University：Natural Sciences,2012,31(1):121-123.(in Chinese)

[10] 李先明,王占军.小井沟水库围堰溃决洪水演进预测分析 [J].水力发电,2013,39(1):20-22.

Li X M,Wang Z J.Prediction on cofferdam break flood routing of Xiaojinggou reservoir [J].Water Power,2013,39(1):20-22.(in Chinese)

[11] 王领元.应用MIKE对河流一、二维的数值模拟 [D].辽宁大连:大连理工大学,2007.

Wang L Y.The numerical simulation of l-D and 2-D river on MIKE [D].Dalian，Liaoning:Dalian University of Technology,2007.(in Chinese)

[12] 李 洋.GIS 支持下的黑河金盆水库溃坝洪水数值模拟 [D].西安:西安理工大学,2007.

Li Y.The Numerical Simulation of dam break floods based on GIS in Jinpen reservoir [D].Xi’an：Xi’an University of Technology,2007.(in Chinese)