一种河道溃决洪水一、二维数学模型耦合分析方法

摘要

本发明提供了一种河道溃决洪水一、二维数学模型耦合分析方法。采集相应数据后，首先采用四边形或三角形非结构网格离散溃口出流淹没区和受溃口出流影响明显的河道段，其余河段采用一维断面数据概化。一维模型采用经典的Presimann隐格式求解，二维模型采用具有激波捕捉能力的Godunov格式求解。为避免数值震荡，保证一、二维模型连接处质量守恒和动量守恒，本方法提出的连接处为一过渡区域，在过渡区域内既进行一维模型计算，也进行二维模型计算。溃口启动后，溃口出流淹没区和河道二维计算区自动成为一体化二维计算区域，溃口出流过程的计算由二维模型完成。这种处理方法能够充分考虑溃口处河道内外的水动力响应过程，提高溃口出流过程的数值模拟精度，克服了传统方法中仅采用堰流公式计算溃口出流过程，不能考虑河道内外动量交换的缺点。

说明书

技术领域

本方法涉及水利工程领域，尤其涉及防洪减灾领域，具体为一种河道溃决洪水一、二维数学模型耦合分析方法。

背景技术

堤防是沿江河、湖泊、海洋的岸边或蓄滞洪区、圩区、水库库区的周边修建的防止洪水漫溢或风暴袭击的挡水建筑物，是防洪工程体系中的重要组成部分。历史上，我国堤防溃决的记录是惊人的，近些年，随着国家对水利工程的投入不断增加，大江大河的堤防防洪标准和工程质量有了明显提高，但是溃堤事件还是时有发生。

目前，在进行溃堤洪水的分析计算时，有两种常用的方法。一种是水文学方法，该方法基于水量平衡的原理，结合地理信息数据，能给出洪水的最终淹没范围，但是该方法一般用于淹没范围相对较小和较封闭的区域，不能提供详细的水动力淹没过程。另一种是水动力学的方法，由于一维模型效率高，所需地形资料相对较少，同时能够很方便地处理一些水工建筑物，被广泛地用于河道洪水模拟中，但在处理一些如无固定路径的地表洪水演进等具有明显二维属性的水流时则具有很大的局限性。二维模型在处理复杂水流上具有明显的优势，精度更高，但是也存在如计算效率较低等缺点。一维模型和二维模型具有各自的优势，将两者进行耦合来模拟溃堤洪水是优化的选择。

目前，一、二维耦合的溃堤数学模型技术已经发展的比较稳定，基本上耦合位置都是设在溃口处，河道完全用一维模型模拟，淹没区用二维模型模拟，这种处理方法有一个很突出的问题，就是当溃决发生后，溃口附近的河道内会出现一个非常复杂的水动力响应过程，在溃口处沿着垂直河道的方向看，临近溃口的河道水面下降迅速，远离溃口的一侧河道水面下降缓慢，尤其是当河道较宽时，溃口对岸的水位甚至不会发生非常明显的变化，这个现象单纯靠一维模型中设定的河道溃口处能量守恒、水量守恒等条件假设来进行求解，会均化溃口附近河道内水动力变化过程，给最终的计算结果带来较大的数值误差，另外，这种耦合方式在退水阶段时常会出现数值震荡问题，为了消除数值震荡，会人为采取一些数值光滑手段，这样也会给计算结果带来较大误差。

发明内容

本发明的目的是提供一种河道溃决洪水一、二维数学模型耦合分析方法，该方法区别于传统的一维河道加二维淹没区的耦合计算方式，将受溃口出流影响明显的上下游附近河道包括溃决淹没区一起采用二维数学模型进行模拟，其它河段采用一维数学模型进行模拟，一、二维模型耦合的界面设置在河道内而不是传统的耦合方式设置在溃口处。该方法能够考虑溃口附近河道内外详细的水动力响应过程，提高溃堤洪水计算的模拟精度。

本发明的目的是通过以下方案实现的：

本方法为一种溃堤洪水一、二维数学模型耦合分析方法，该方法针对现有的一、二维堤防溃决洪水数学模型存在的问题，提出将一、二维数学模型的耦合界面设置在河道内(溃口出流影响明显的河段上下游位置处)而不是传统的溃口位置处，能够克服传统一、二维耦合模型对溃口处河道水力响应过于简化的缺点，提升整个溃堤洪水过程的计算精度。该方法需要采集河道地形数据和溃堤淹没区地形数据，溃口信息数据以及河道上、下游洪水信息数据，具体包括以下步骤：

1)设置溃口位置，确定溃口出流影响明显的河道范围。溃口出流淹没区和受溃口出流影响明显的河道段采用二维四边形或三角形非结构网格进行离散，河道其余部分采用一维断面数据概化离散。标记溃口所在位置处的二维网格单元，其网格的溃决属性赋值为“1”，其余网格溃决属性赋值为“0”。

2)明确一维数学模型和二维数学模型耦合时的连接断面和连接网格。为避免在一、二维耦合连接处出现数值震荡，在连接处设置过渡区域，该过渡区域一般可选择一、二维模型连接处2个临近断面之间的区域。为便于表述，一维河道溃口上游段末端两个断面的编号记为I-1和I，一维河道溃口下游段起始端两个断面的编号记为J-1和J，在断面I-1和I之间以及断面J-1和J之间的区域为过渡区域，过渡区域内同样要进行网格剖分。

3)模型初始化。给一维河道断面和二维网格单元设置糙率和水力条件初始值，二维网格高程赋值。

4)确定计算的时间步长dt。一维河道模型采用Presimann隐格式有限差分法求解，对时间步长要求宽松；二维模型采用Godunov显格式有限体积法求解，时间步长受CFL(Courant-Friedrichs-Lewy)条件限制。为保证一、二维模型交互同步性和整体稳定性，计算时间步长统一取为二维模型的计算步长。时间步长确定后，获取t+dt时刻河道首末断面的边界条件信息。

5)二维模型计算。在溃口上游段，一维模型将I-1连接断面t时刻的水位和流量值传递给与其相连接的二维连接单元，二维连接单元以此流量和水位值作为入流边界进行计算；在溃口下游段，一维模型将J连接断面t时刻的水位值和流量值传递给与其相连接的二维连接单元，二维连接单元以此流量和水位值作为出流边界进行计算。二维模型计算一个时间步长dt结束后，所有二维网格单元的水力要素值更新到t+dt时刻。

6)一维模型计算。返回与I断面相连的二维连接单元t+dt时刻的水力要素值给连接断面I，作为溃口上游段一维模型的出流边界进行计算；返回与J-1断面相连的二维连接单元t+dt时刻的水力要素值给连接断面J-1，作为溃口下游段一维模型的入流边界进行计算。一维模型计算一个时间步长dt后，得到t+dt时刻一维河道各断面的水位和流量值。

7)是否启动溃口。不失一般性，假定溃口为瞬时溃决模式，当溃口处达到溃决条件时，将溃口所在位置处溃口属性为“1”的二维网格单元高程改变至溃口真实的溃决高程，溃决淹没计算自动开始，溃口处的出流过程由二维模型自动计算。溃口一旦启动，溃口出流淹没区和河道二维计算区自动成为一体化二维计算区域。

8)更新t＝t+dt，重复步骤4)～(7)直至计算结束。

进一步的，步骤1)中，在确定溃口出流影响明显的河道范围时需统筹考虑河道的地形特征、河道水流特征、溃口宽度等要素。在初始建模阶段，该范围取值不应低于溃口宽度的20倍，即溃口上游受溃口出流影响明显的河道范围至少为溃口宽度的10倍，溃口下游受溃口出流影响明显的河道范围也至少为溃口宽度的10倍。根据模型的初步计算结果，可适当调整该范围的大小，调整依据为：该河道影响范围两端的二维水流状态并不明显时，则说明所确定的范围大小是合适的。

进一步的，步骤2)中，一维、二维数学模型连接处不是传统方法上用的一个界面(河道和溃口连接位置)而是河道内部远离溃口的两个过渡区域，在I-1和I断面之间以及J-1和J断面之间既要进行一维计算，也要进行二维计算。该处理方法可保证整体水量、动量守恒，避免模型连接处出现数值震荡。

进一步的，步骤4)确定计算的时间步长dt统一取为二维模型的计算步长，一维河道模型采用Presimann隐格式有限差分法求解，控制方程组为扩展的Saint-Venant方程组，如式(1)、(2)所示：

质量方程

动量方程

式中：X为河长、t为时间；A为过水面积；Q为断面流量；Z为断面水位；α为动量修正系数；K为流量模数；q_l为旁侧入流；v_x为入流沿水流方向的速度分量，若旁侧入流垂直于干流，则v_x＝0；

二维模型采用Godunov显格式有限体积法求解，采用的控制方程组为守恒型二维浅水方程组，如式(3)所示：

其中：

h为水深，t为时间，u，v分别为x，y方向的流速；

分别为x，y方向的底坡项，其中Z_b为底高程；

分别为x，y方向的摩阻坡降，其中n为曼宁糙率系数。

进一步的，步骤6)中，二维模型在给一维模型传递水力要素值时，需要将二维模型的水力要素值u，v和h转为一维模型所需的过水断面流量值和水位值。转化方法如下：

二维模型采用的浅水方程的连续方程如下式所示：

式中：h为水深；u，v分别为x，y方向的流速；t为时间。

当采用Godunov显格式有限体积法进行求解时，连续方程的数值通量值即为通过单元边的法向单宽流量值，因此将I和J-1断面处的二维单元边乘以该单元边处的流量数值通量，然后求和即可得到断面I和断面J-1处的断面流量值，严格满足模型交互处的质量守恒。将与I和J-1断面相连接的二维单元水深值加上各自的底高程，然后取平均值即为断面I和断面J-1处的断面水位值。

进一步的，步骤7)中溃口一旦启动，溃口出流淹没区和河道二维计算区自动成为一体化二维计算区域，溃口出流过程的计算由二维模型完成，这种处理方法能够充分考虑溃口处河道内外的水动力响应过程，提高溃口出流过程的数值模拟精度，克服了传统方法中仅采用堰流公式计算溃口出流过程，不能考虑河道内外动量交换的缺点。

本发明的有益效果：

在进行溃堤水流一、二维数学模型耦合分析计算时，将一、二维模型的耦合位置设置在河道内(溃口出流影响明显的河段上下游位置处)而不是传统的溃口位置处，传统的一、二维数学模型耦合位置为一界面，本方法提出的耦合位置为一过渡区域，这种处理方式能够有效保证整体计算的质量守恒和动量守恒，避免数值震荡现象的发生，另外，本方法可充分考虑溃口附近河道内外复杂的水力条件变化，使得模拟结果更接近真实情景。综上，本方法能够有效提高现有一、二维耦合溃堤数学模型的计算精度。

附图说明

图1为本发明的溃堤洪水一、二维数学模型耦合分析方法流程图；

图2为本发明的溃堤洪水一、二维数学模型耦合分析方法示意图。

具体实施方式

下面结合附图1和附图2对本发明作进一步详细说明。

本发明提供的是一种河道溃决洪水一、二维数学模型耦合分析方法。该方法中，一维、二维数学模型连接处不是传统方法上用的一个界面(河道和溃口连接处)而是河道内部的两个过渡区域，在该过渡区域内既要进行一维计算，也要进行二维计算，可有效避免数值震荡现象出现，充分保证数值计算的质量和动量守恒。溃口一旦启动，溃口出流淹没区和河道二维计算区自动成为一体化二维计算区域，溃口出流过程自动由二维模型计算，该种处理方法可充分考虑河道内外水力条件的变化，提高模拟精度。该方法主要包括以下具体步骤：

(1)设置溃口位置，确定溃口出流影响明显的河道范围。溃口出流淹没区和受溃口出流影响明显的河道段采用二维四边形或三角形非结构网格进行离散，河道其余部分采用一维断面数据概化离散。标记溃口所在位置处的二维网格单元，其网格的溃决属性赋值为“1”，其余网格溃决属性赋值为“0”。在溃口处，为更细致的描绘溃口形状，可适当对网格进行加密剖分处理。

在确定溃口出流影响明显的河道范围时需统筹考虑河道的地形特征、河道水流特征、溃口宽度等要素。在初始建模阶段，该范围取值不应低于溃口宽度的20倍，即溃口上游受溃口出流影响明显的河道范围至少为溃口宽度的10倍，溃口下游受溃口出流影响明显的河道范围也至少为溃口宽度的10倍。根据模型的初步计算结果，可适当调整该范围的大小，调整依据为：该河道影响范围两端的二维水流状态并不明显时，则说明所确定的范围大小是合适的。

(2)明确一维数学模型和二维数学模型耦合时的连接断面和连接网格。为避免在一、二维耦合连接处出现数值震荡，在连接处设置过渡区域，该过渡区域一般可选择一、二维模型连接处2个临近断面之间的区域。为便于表述，一维河道溃口上游段末端两个断面的编号记为I-1和I，一维河道溃口下游段起始端两个断面的编号记为J-1和J，在断面I-1和I之间以及断面J-1和J之间的区域为过渡区域，过渡区域内同样要进行网格剖分。

(3)模型初始化。给一维河道断面和二维网格单元设置糙率和水力条件初始值，二维网格高程赋值。

(4)确定计算的时间步长dt。一维河道模型采用Presimann隐格式有限差分法求解，控制方程组为扩展的Saint-Venant方程组，如下式所示：

质量方程

动量方程

式中：x、t为河长、时间；A为过水面积；Q为断面流量；Z为断面水位；α为动量修正系数；K为流量模数；q_l为旁侧入流；v_x为入流沿水流方向的速度分量，若旁侧入流垂直于干流，则v_x＝0。

采用Preissmann四点隐式差分格式对上述方程组进行离散后得到一个离散化的非线性方程组，采用具有较高收敛速度的Newton-Raphson迭代法来直接求解非线性代数方程组，详细过程可参见以下著作(张大伟，基于Godunov格式的地方溃决水流数值模拟，2014，中国水利水电出版社)。

二维模型采用Godunov显格式有限体积法求解，采用的控制方程组为守恒型二维浅水方程组，如下式所示：

其中：

h为水深，t为时间，u，v分别为x，y方向的流速。

分别为x，y方向的底坡项，其中Z_b为底高程。

分别为x，y方向的摩阻坡降，其中n为曼宁糙率系数。

详细的二维模型构建过程可以参见如下著作(张大伟，基于Godunov格式的地方溃决水流数值模拟，2014，中国水利水电出版社)。

由于采用隐格式进行求解，所以一维模型对时间步长的要求比较宽松。二维模型为显格式求解，时间步长选取受CFL(Courant-Friedrichs-Lewy)条件限制。为保证一、二维模型交互同步性和整体稳定性，计算时间步长统一取为二维模型的计算步长。时间步长确定后，获取t+dt时刻河道首末断面的边界条件信息。

(5)二维模型计算。在溃口上游段，一维模型将I-1连接断面t时刻的水位和流量值传递给与其相连接的二维连接单元，二维连接单元以此流量和水位值作为入流边界进行计算；

在溃口下游段，一维模型将J连接断面t时刻的水位值和流量值传递给与其相连接的二维连接单元，二维连接单元以此流量和水位值作为出流边界进行计算。二维模型计算一个时间步长dt结束后，所有二维网格单元的水力要素值更新到t+dt时刻。

(6)一维模型计算。返回与I断面相连的二维连接单元t+dt时刻的水力要素值给连接断面I，作为溃口上游段一维模型的出流边界进行计算；返回与J-1断面相连的二维连接单元t+dt时刻的水力要素值给连接断面J-1，作为溃口下游段一维模型的入流边界进行计算。

二维模型在给一维模型传递水力要素值时，需要将二维模型的水力要素值u，v和h转为一维模型所需的过水断面流量值和水位值。转化方法如下：

二维模型采用的浅水方程的连续方程如下式所示：

式中：h为水深；u，v分别为x，y方向的流速，t为时间。

二维模型的水力要素值为网格中心点处的水深和流速，一维模型需要的边界为过水断面处的水位和流量。所以将二维模型计算得到的水力要素值提供给一维模型时，需要进行以下转换：当采用Godunov显格式有限体积法进行求解时，连续方程的数值通量值即为通过单元边的法向单宽流量值，因此将I和J-1断面处的二维单元边乘以该单元边处的流量数值通量，然后求和即可得到断面I和断面J-1处的断面流量值，严格满足模型交互处的质量守恒。将与I和J-1断面相连接的二维单元水深值加上各自的底高程，然后取平均值即为断面I和断面J-1处的断面水位值。一维模型计算一个时间步长dt后，得到t+dt时刻一维河道各断面的水位和流量值。

(7)是否启动溃口。不失一般性，假定溃口为瞬时溃决模式，当溃口处达到溃决条件时，将溃口所在位置处溃口属性为“1”的二维网格单元高程改变至溃口真实的溃决高程，溃决淹没计算自动开始，溃口处的出流过程由二维模型自动计算。溃口一旦启动，溃口出流淹没区和河道二维计算区自动成为一体化二维计算区域。溃口出流过程的计算由二维模型完成，这种处理方法能够充分考虑溃口处河道内外的水动力响应过程，提高溃口出流过程的数值模拟精度，克服了传统方法中仅采用堰流公式计算溃口出流过程，不能考虑河道内外动量交换的缺点。

(8)更新t＝t+dt，重复步骤(4)～(7)直至计算结束。

上述的实施例仅是本发明的部分体现，并不能涵盖本发明的全部，在上述实施例以及附图的基础上，本领域技术人员在不付出创造性劳动的前提下可获得更多的实施方式，因此这些不付出创造性劳动的前提下获得的实施方式均应包含在本发明的保护范围内。