本发明公开了一种包含库区的土石坝漫顶溃决洪水数值模拟方法。采用三角形非结构网格对库区和下游淹没区进行剖分，坝体所在区域不参与剖分。将溃口出流量作为库区的出流边界和下游淹没区的入流边界，分别进行库区和下游淹没区的水动力演进计算，出流边界和入流边界所对应的坝前单元和坝后单元根据溃口的实时宽度动态调整。该方法不仅能够提供下游的淹没形态，而且能够提供库区内的水流运动过程，从而扩大溃坝模型的适用范围。在计算溃口出流时，使用跟溃口直接相连的坝前单元的平均水位，而不是整个库区的平均水位，可显著提高溃口出流过程的准确性，从而完善现有的土石坝溃决洪水计算方法，提高土石坝溃决洪水计算的数值精度。
1.一种包含库区的土石坝漫顶溃决洪水数值模拟方法，其特征在于：将库区包含在整个计算区域内，土石坝溃口出流过程的模拟与库区和下游淹没区洪水演进模拟统一进行，具体步骤如下：
1)获取所需的基础数据：所述基础数据包括库区水下地形数据、坝体相关数据、溃坝淹没区地形数据、整个计算区域DOM地貌数据、入库洪水过程数据；所述坝体相关数据包括坝高、坝长、坝顶宽度、迎水面和背水面坡度、坝体材料粘性系数、粒径级配；
2)计算区域网格离散：库区与坝体相连接的边定义为坝前边，下游淹没区与坝体相连接的边定义为坝后边，坝前边与坝后边长度相等，相互对应；坝前边和坝后边作为网格剖分的控制边，坝前边与坝后边的离散方式相同，离散尺度控制在5m；采用三角形非结构网格对库区和下游淹没区分别离散，坝体所在的区域不参与网格离散，库区和下游淹没区网格单元统一编码；离散后，坝前边对应的网格单元称为坝前单元，坝后边对应的网格单元称为坝后单元，坝前单元和坝后单元数量相同，一一对应；
3)计算区域赋值：采用库区和下游淹没区地形数据对网格单元型心位置进行高程插值；根据地貌数据对网格单元糙率赋值；对网格单元的水力变量赋初值；将入库洪水过程设置为库区的入流边界；下游淹没区的出流边界设置为自由出流；
4)溃口设置：设定溃口形状为矩形，溃口初始宽度设置为0.5m，初始深度设置为0.2m；
不失一般性，溃口初始位置设置在坝前边的中间位置，此时溃口对应的坝前单元编号为U0，以U0单元为中心，顺水流方向左侧坝前单元编号依次为U1，U3，U5……Un，其中n为奇数，顺水流方向右侧坝前单元编号依次为U2，U4，U6……Um，其中m为偶数；与溃口初始位置对应的坝后单元编号为D0，以D0单元为中心，顺水流左侧坝后单元编号依次为D1，D3，D5……Dn，其中n为奇数，顺水流右侧坝后单元编号依次为D2，D4，D6……Dm，其中m为偶数；按照各坝前单元及坝后单元的编号顺序，与各坝前单元及坝后单元分别对应的坝前边及坝后边的长度依次定义为B0，B1，B2，B3，B4，B5，B6……Bi，i为坝前边或坝后边离散后线段的总条数，i为自然数；
5)确定时间步长dt：由各网格单元的水力变量值，获取满足计算稳定CFL条件的时间步长dt；
6)溃口出流过程计算：根据溃口出流公式计算t时刻溃口的出流量，溃口出流公式如式(1)所示：Qb＝kbB(Zu-Zd)1.5 (1)
其中：Qb为溃口出流量，kb为流量系数，B为溃口宽度，Zu为与溃口直接相连的所有坝前单元的平均水位值，Zd为溃口底部高程；
根据de Vries输沙率公式来计算溃口处的泥沙输移量，以确定下一时刻t+dt的溃口宽度和深度，计算溃口发展；
7)洪水演进过程计算：洪水演进计算包括库区洪水演进计算和下游淹没区洪水演进计算，采用守恒型的二维浅水方程组描述洪水运动；使用溃口流量Qb驱动库区和下游淹没区的洪水演进计算，在库区洪水演进计算中Qb作为出流边界处理，在下游淹没区洪水演进计算中Qb作为入流边界处理；
根据溃口在当前时刻的宽度B来确定出流边界所对应的坝前单元和入流边界对应的坝后单元：在选择坝前单元时，如果B小于B0，出流边界所对应的坝前单元为U0；当B大于B0时，以U0为中心，每一次在原来已选择的坝前单元的基础上再向两侧各选取1个坝前单元，直到所选的所有坝前单元对应的坝前边之和大于等于溃口宽度B时结束，选中的坝前单元作为当前时刻出流边界所在单元参与本时段溃口流量Qb的分配；与选中的坝前单元对应的坝后单元作为当前时刻入流边界所在单元参与本时段溃口流量Qb的分配；
8)更新t＝t+dt，重复步骤5)～7)直至计算结束。
2.根据权利要求1所述的包含库区的土石坝漫顶溃决洪水数值模拟方法，其特征在于：
其中，dZd为时段dt内溃口的冲刷深度，cs和βs为输沙率公式系数，Qb为溃口出流量，Ab为溃口的过水断面面积；设定溃口的冲刷深度与侧向展宽的宽度相同，溃口向两侧对称展宽。
3.根据权利要求1所述的包含库区的土石坝漫顶溃决洪水数值模拟方法，其特征在于：
步骤7)中，描述洪水运动所采用的方程组为守恒型二维浅水方程组，采用具有良好激波捕捉能力的Godunov格式离散该方程组，采用的守恒型二维浅水方程组如下式所示：式中：
h为水深，u，v分别为x，y方向的流速，t为时间， 分别x，y方向的
坡度，Zb为地面高程，g为重力加速度， 分别为x，y方向的
一种包含库区的土石坝漫顶溃决洪水数值模拟方法
技术领域
[0001] 本发明涉及水利工程领域，尤其涉及防洪减灾领域，具体为一种包含库区的土石坝漫顶溃决洪水数值模拟方法。
背景技术
[0002] 我国是世界上修建水库大坝最多的国家，根据第一次全国水利普查工作统计，我国现有水库大坝98002座，其中土石坝作为应用最广泛的坝型，约占已建坝总数的95％以上。已经修建的土石坝中，绝大多数为小型水库，由于特殊的历史原因，这些小型土石坝在修建时多缺乏严格的设计和施工程序，当遭遇超标准洪水时，存在比较较大的溃决风险。水库大坝一旦溃决失事，将会给下游带来巨大的灾难。在土石坝的溃决失事方式中，最常见的为洪水漫顶溃决，洪水一旦漫顶，会迅速在坝后形成冲刷，直至溃口稳定。另外值得说明的是，土石坝溃决洪峰流量远大于水库的正常泄洪流量，所以在溃坝分析计算时，一般可忽略溢洪道泄洪或底孔泄洪的影响。
[0003] 对土石坝溃决洪水分析计算来讲，非常重要的一步是计算溃口的出流过。溃口出流计算的方法总体来说可以分为两大类：一类为采用溃坝的历史数据资料得出的回归方程来计算出流过程，这类回归过程一般是建立溃口出流过程与库容和坝高的相关关系的表达式。另一类为采用带有溃口发展物理机制的理论或半理论方法来预测出流过程，在这其中又可分为两小类，一类为参数型溃口模型，另一类为物理-参数型溃口模型。所谓参数型溃口模型是指通过设定溃口尺寸，形状和发展时间，利用堰流公式来计算溃坝洪水出流过程，这些模型通常采用相对简单的参数化经验方法来估算溃口发展过程。所谓物理-参数型溃口模型指的是溃口的发展考虑了溃口冲刷的物理机制，通过泥沙输移公式来计算溃口侵蚀过程。这类模型相对较多，代表了土石坝溃口出流模型的一个发展方向。
[0004] 现在的土石坝溃决计算实际上人为分成了两个步骤，首先是溃口出流过程的计算，得到完整的溃口出流过程后，再进行下一步洪水演进计算。然而，真实的溃决过程是库区水位降低、溃口发展和下游洪水演进是同时进行的，也是相互影响的，溃口前局部水位的高低直接影响着溃口的发展，溃口的发展决定了溃口出流量的大小，出流量的大小决定了下游淹没形势的不同，因此，将库区包含在模型内，进行库区、溃口和下游淹没区一体化的模拟是非常有必要的，能够提高当前土石坝溃决洪水的模拟精度。
发明内容
[0005] 本发明的目的在于提供一种包含库区的土石坝漫顶溃决数值模拟方法，在每个时间步长内同时完成库区洪水演进计算、溃口出流计算和下游淹没过程的计算，充分考虑土石坝溃决过程的整体性，提高现有土石坝溃决洪水分析计算的精度。
[0006] 本发明是通过以下技术方案实现的：
[0007] 将库区包含在整个计算区域内，土石坝溃口出流过程的模拟与库区和下游淹没区洪水演进模拟统一进行，具体步骤如下：
[0008] 1)获取所需的基础数据：所述基础数据包括库区水下地形数据(比例尺1:5000及以上)；坝体相关数据，包含坝高、坝长、坝顶宽度、迎水面和背水面坡度、坝体材料粘性系数、粒径级配等；溃坝淹没区地形数据(比例尺1:10000及以上)；获取整个计算区域DOM地貌数据(分辨率2m及以上)；入库洪水过程数据。
[0009] 2)计算区域网格离散：库区与坝体相连接的边定义为坝前边，下游淹没区与坝体相连接的边称为坝后边，坝前边与坝后边长度相等，相互对应；坝前边和坝后边作为网格剖分的控制边，离散方式相同，离散尺度控制在5m左右；采用三角形非结构网格对库区和下游淹没区分别离散，坝体所在的区域不参与网格离散，库区和下游淹没区网格单元统一编码；离散后，坝前边对应的网格单元称为坝前单元，坝后边对应的网格单元为坝后单元，坝前单元和坝后单元数量相同，一一对应。
[0010] 3)计算区域赋值：采用库区和下游淹没区地形数据对网格单元型心位置进行高程插值；根据地貌数据对网格单元糙率赋值；对网格单元的水力变量赋初值；将入库洪水过程设置为库区的入流边界；下游淹没区的出流边界设置为自由出流。
[0011] 4)溃口设置：设定溃口形状为矩形，溃口初始宽度设置为0.5m，初始深度设置为0.2m；不失一般性，溃口初始位置设置在坝前边的中间位置，此时溃口对应的坝前单元编号为U0，以U0单元为中心，顺水流方向左侧坝前单元编号依次为U1，U3，U5……Un，其中n为奇数，，顺水流方向右侧坝前单元编号依次为U2，U4，U6……Um，其中m为偶数；同样的，与溃口初始位置对应的坝后单元编号为D0，以D0单元为中心，顺水流左侧坝后单元编号依次为D1，D3，D5……Dn，其中n为奇数，顺水流右侧坝后单元编号依次为D2，D4，D6……Dm，其中m为偶数；按照各坝前单元(坝后单元)的编号顺序，与各坝前单元(坝后单元)对应的坝前边(坝后边)的长度依次定义为B0，B1，B2，B3，B4，B5，B6……Bi，i为坝前边或坝后边离散后线段的总条数，i为自然数。
[0012] 5)确定时间步长dt：由各网格单元的水力变量值，获取满足计算稳定CFL条件的时间步长dt。
[0013] 6)溃口出流过程计算：根据溃口出流公式计算t时刻溃口的出流量，公式如下式所示：
[0014] Qb＝kbB(Zu-Zd)1.5 (1)
[0015] 其中：Qb为溃口出流量，Kb为流量系数，B为溃口宽度，Zu为与溃口直接相连的所有坝前单元的平均水位值，Zd为溃口底部高程。
[0016] 根据de Vries输沙率公式来计算溃口处的泥沙输移量，以确定下一时刻t+dt的溃口宽度和深度，计算溃口发展。
[0017] 7)洪水演进过程计算：洪水演进计算包括库区洪水演进计算和下游淹没区洪水演进计算，采用守恒型的二维浅水方程组描述洪水运动；使用溃口流量Qb驱动库区和下游淹没区的洪水演进计算，在库区洪水演进计算中Qb作为出流边界处理，在下游淹没区洪水演进计算中Qb作为入流边界处理；
[0018] 根据溃口在当前时刻的宽度B来确定出流边界所对应的坝前单元和入流边界对应的坝后单元。在选择坝前单元时，，如果B小于B0，出流边界所对应的坝前单元为U0；当B大于B0时，以U0为中心，每一次在原来已选择的坝前单元的基础上再向两侧各选取1个坝前单元，直到所选的所有坝前单元对应的坝前边之和大于等于溃口宽度B时结束，这些选中的坝前单元作为当前时刻出流边界所在单元参与本时段溃口流量Qb的分配；同理，与这些选中的坝前单元对应的坝后单元作为当前时刻入流边界所在单元参与本时段溃口流量Qb的分配。
[0019] 8)更新t＝t+dt，重复步骤5)～7)直至计算结束。
[0020] 进一步的，步骤6)中溃口发展的计算公式如下：
[0021]
[0022] 其中，dZd为时段dt内溃口的冲刷深度，cs和βs为输沙率公式系数，Qb为溃口出流量，Ab为溃口的过水断面面积。设定溃口的冲刷深度与侧向展宽的宽度相同，溃口向两侧对称展宽。
[0023] 进一步的，步骤7)中，描述洪水运动所采用的方程组为守恒型二维浅水方程组，采用具有良好激波捕捉能力的Godunov格式离散该方程组。采用的守恒型二维浅水方程组如下式所示：
[0024]
[0025] 式中：
[0026] h为水深，u，v分别为x，y方向的流速，t为时间， 分别x，y方向的坡度，Zb为地面高程，g为重力加速度， 分别
为x，y方向的摩阻项，其中n为Manning糙率系数。
[0027] 本发明的有益效果：
[0028] 在进行土石坝漫顶溃决模拟时，通过将溃口出流过程分别作为上下游区域的出流和入流边界条件的方法将库区、溃口和下游淹没区统一起来建立一体化的土石坝溃决模拟模型，该方法不仅能够提供下游的淹没形态，而且能够提供库区内的水流运动过程，在计算溃口出流时，使用溃口前的局部水位，而不是整个库区的平均水位，可显著提高溃口出流过程的准确性，从而整体上完善现有的土石坝溃决洪水计算方法，提高土石坝溃决洪水计算的数值精度。
附图说明
[0029] 图1为本发明的一种包含库区的土石坝漫顶溃决洪水数值模拟方法流程图；
[0030] 图2为本发明的一种包含库区的土石坝漫顶溃决洪水数值模拟方法示意图。
具体实施方式
[0031] 下面结合附图1和附图2对本发明作进一步说明。
[0032] 本发明为一种包含库区的土石坝漫顶溃决洪水数值模拟方法，将库区包含在整个计算区域内，土石坝溃口出流过程的模拟与库区和下游淹没区洪水演进模拟统一进行，具体步骤如下：
[0033] 1)获取所需的基础数据：所述基础数据包括库区水下地形数据(比例尺1:5000及以上)；坝体相关数据，包含坝高、坝长、坝顶宽度、迎水面和背水面坡度、坝体材料粘性系数、粒径级配等；溃坝淹没区地形数据(比例尺1:10000及以上)；获取整个计算区域DOM地貌数据(分辨率2m及以上)，入库洪水过程数据。
[0034] 2)计算区域网格离散：库区与坝体相连接的边定义为坝前边，下游淹没区与坝体相连接的边称为坝后边，坝前边与坝后边长度相等，相互对应；坝前边和坝后边作为网格剖分的控制边，离散方式相同，离散尺度控制在5m左右；采用三角形非结构网格对库区和下游淹没区分别离散，坝体所在的区域不参与网格离散，库区和下游淹没区网格单元统一编码；离散后，坝前边对应的网格单元称为坝前单元，坝后边对应的网格单元为坝后单元，坝前单元和坝后单元数量相同，一一对应。
[0035] 3)计算区域赋值：采用库区和下游淹没区地形数据对网格单元型心位置进行高程插值；根据地貌数据对网格单元糙率赋值；对网格单元的水力变量赋初值；将入库洪水过程设置为库区的入流边界；下游淹没区的出流边界设置为自由出流。
[0036] 4)溃口设置：设定溃口形状为矩形，溃口初始宽度设置为0.5m，初始深度设置为0.2m；不失一般性，溃口初始位置设置在坝前边的中间位置，此时溃口对应的坝前单元编号为U0，以U0单元为中心，顺水流方向左侧坝前单元编号依次为U1，U3，U5……Un，其中n为奇数，顺水流方向右侧坝前单元编号依次为U2，U4，U6……Um，其中m为偶数；同样的，与溃口初始位置对应的坝后单元编号为D0，以D0单元为中心，顺水流左侧坝后单元编号依次为D1，D3，D5……Dn，其中n为奇数，，顺水流右侧坝后单元编号依次为D2，D4，D6……Dm，其中m为偶数；按照各坝前单元(坝后单元)的编号顺序，与各坝前单元(坝后单元)对应的坝前边(坝后边)的长度依次定义为B0，B1，B2，B3，B4，B5，B6……Bi，i为坝前边或坝后边离散后线段的总条数，i为自然数。
[0037] 5)确定时间步长dt：由各网格单元的水力变量值，获取满足计算稳定CFL条件的时间步长dt。
[0038] 6)溃口出流过程计算：根据溃口出流公式计算t时刻溃口的出流量，公式如下式所示：
[0039] Qb＝kbB(Zu-Zd)1.5 (1)
[0040] 其中：Qb为溃口出流量，Kb为流量系数，B为溃口宽度，Zu为与溃口直接相连的所有坝前单元的平均水位值，以图2为例，此时的Zu为U0，U1和U2三个坝前单元水位的平均值，Zd为溃口底部高程。
[0041] 根据de Vries输沙率公式来计算溃口处的泥沙输移量，以确定下一时刻t+dt的溃口宽度和深度，溃口发展的计算公式如下：
[0042]
[0043] 其中，dZd为时段dt内溃口的冲刷深度，cs和βs为输沙率公式系数，Qb为溃口出流量，Ab为溃口的过水断面面积。设定溃口的冲刷深度与侧向展宽的宽度相同，溃口向两侧对称展宽。当达到设定的溃口展宽宽度后，溃口发展达到稳定状态，溃口的形状不在发生变化。
[0044] 7)洪水演进过程计算：洪水演进过程计算包括库区洪水演进计算和下游淹没区洪水演进计算，采用守恒型的二维浅水方程组描述洪水运动。守恒型二维浅水方程组如下式所示：
[0045]
[0046] 式中： h为水深，u，v分别为x，y方向的流速，t为时间， 分别x，y方向的坡度，Zb为
地面高程，g为重力加速度， 分别为x，y方向的摩
阻项，其中n为Manning糙率系数。
[0047] 对以上浅水方程组进行非结构离散，采用Roe格式的近似Riemann解计算通过非结构单元界面的数值通量，对底坡项采用特征分解的方法处理，摩阻项采用半隐的方法处理。详细的数值求解过程可参见如下著作P85-P90页(张大伟.基于Godunov格式的堤坝溃决水流数值模拟,2014,中国水利水电出版社)。
[0048] 使用溃口流量Qb驱动库区和下游淹没区的洪水演进计算，在库区洪水演进中Qb作为出流边界处理，在下游淹没区洪水演进计算中Qb作为入流边界处理。根据溃口在当前时刻的宽度B来确定出流边界所对应的坝前单元和入流边界对应的坝后单元。在选择坝前单元时，如果B小于B0，出流边界所对应的坝前单元为U0；当B大于B0时，以U0为中心，每一次在原来已选择的坝前单元的基础上再向两侧各选取1个坝前单元，直到所选的所有坝前单元对应的坝前边之和大于等于溃口宽度B时结束，如图2所示，此时的溃口宽度B明显大于B0，则以U0为中心，选取U1和U2后，B0+B1+B2之和已大于B，选取过程结束，选中的坝前单元为U0，U1和U2，这些选中的坝前单元作为当前时刻出流边界所在单元参与本时段溃口流量Qb的分配；同理，与这些选中的坝前单元对应的坝后单元作为当前时刻入流边界所在单元参与本时段溃口流量Qb的分配。详细的Qb分配处理方法可参见如下著作P92页(张大伟.基于Godunov格式的堤坝溃决水流数值模拟,2014,中国水利水电出版社)。
[0049] 8)更新t＝t+dt，重复步骤5)～7)直至计算结束。
[0050] 上述的实施例仅是本发明的部分体现，并不能涵盖本发明的全部，在上述实施例以及附图的基础上，本领域技术人员在不付出创造性劳动的前提下可获得更多的实施方式，因此这些不付出创造性劳动的前提下获得的实施方式均应包含在本发明的保护范围内。
