CN108021780B - 一种基于无规则非结构网格模型的山洪动态仿真方法 - Google Patents

一种基于无规则非结构网格模型的山洪动态仿真方法 Download PDF

Info

Publication number
CN108021780B
CN108021780B CN201810067029.XA CN201810067029A CN108021780B CN 108021780 B CN108021780 B CN 108021780B CN 201810067029 A CN201810067029 A CN 201810067029A CN 108021780 B CN108021780 B CN 108021780B
Authority
CN
China
Prior art keywords
flood
grid
points
flux
boundary
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810067029.XA
Other languages
English (en)
Other versions
CN108021780A (zh
Inventor
马军建
徐青
熊光亚
郑建兵
余有胜
夏洲
芮钧
胡鸣东
季卫东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nari Technology Co Ltd
Original Assignee
Nari Technology Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nari Technology Co Ltd filed Critical Nari Technology Co Ltd
Priority to CN201810067029.XA priority Critical patent/CN108021780B/zh
Publication of CN108021780A publication Critical patent/CN108021780A/zh
Application granted granted Critical
Publication of CN108021780B publication Critical patent/CN108021780B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于无规则非结构网格模型的山洪动态仿真方法，首先选定模拟区域，收集资料；采用无规则非结构网格针对复杂计算域和不规则地形进行网格剖分；采用二维非恒定流浅水方程组描述水流流动，用有限体积法对和黎曼近似解求解方程组；计算每个网格单元的洪水到达时间、洪水淹没时间、最大洪水淹没水深、最大洪水水位和最大洪水流速；得到洪水到达时间等值图、洪水淹没时间等值图、最大洪水淹没水深等值图，最大的洪水水位等值图和最大洪水流速等值图。本发明实现了复杂计算域和不规则地形的山洪易发山丘区的洪水淹没过程的动态模拟和仿真，实时模拟和计算山洪淹没范围、淹没水深、流场等洪水要素，实现洪水淹没过程的可视化仿真。

Description

一种基于无规则非结构网格模型的山洪动态仿真方法
技术领域
本发明涉及一种山洪动态仿真计算方法，具体涉及一种基于无规则非结构网格模型的山洪动态仿真方法。
背景技术
山洪灾害是山丘区在一定强度或持续降雨下，因特殊的地形地质条件而发生的自然灾害，它具有突发、易发、多发、破坏性大和防治困难的鲜明特点，往往对局部地区造成毁灭性灾害，造成重大的生命财产损失。为了对山洪产生到发展趋势进行动态分析和及时预警，以及预测和分析评估山洪造成的灾害情况，通过基于物理过程的山洪动态演进仿真可有效实现这一目标，并能提供沿程水位、流速、淹没范围等重要信息，实现洪水演进的准确模拟。
目前，在山洪演进研究手段上主要包括水动力模拟、洪水波理论分析、水文统计等。以往河道汇流往往采用传统水文学方法或采用单位线近似，由于忽略了力学因子使得传统方法难以准确描述水流运动的问题，且无法描述洪水溃坝、溃堤后的二维水流演进。水动力模拟能从洪水运动的物理机理出发，利用流体控制微分方程和数值离散算法，准确的模拟洪水波传播运动规律，显示出强大的优越性，已成为防洪减灾工作中的重要手段。
常用的洪水演进水动力模拟方法包括有限体积法、有限差分法、有限元法等，其中有限体积法具有良好的守恒性和计算精度，逐渐成为模拟的主流方法之一。在山洪演进模拟过程中，往往发生流量、水位在极短时间内急剧增加，洪水波波形传播过程中发生间断，这往往要采取特殊的方法做特殊的处理，但在计算格式上仍缺乏和谐性和稳定性，模拟结果不够理想。此外，山洪易发区往往分布在山丘区，属于强无规则型地形，计算过程前处理中网格剖分具有非结构特点，常规的水动力模拟模型难以适应。
发明内容
本发明的目的在于克服现有技术中的不足，针对复杂计算域和强不规则地形的山洪演进模拟问题，本发明提供了一种基于无规则非结构网格模型的山洪动态仿真方法，对强无规则型的地形进行网格加密剖分处理，采用二维非恒定流浅水方程组描述水流流动，应用有限体积法及黎曼近似解对耦合方程组进行数值求解，从而模拟出山洪动态演进过程。
为解决上述技术问题，本发明提供了一种基于无规则非结构网格模型的山洪动态仿真方法，包括以下步骤：
(1)选定模拟区域，山洪易发的典型区域，或易发生垮堤的土堤河段，输入资料数据包括：山洪淹没区域的气候、水文、地形、地貌、地质及水工建筑物基本概况资料，山洪淹没区域的水上数字高程模型DEM地形数据和河道地形数据，山洪淹没区域的分类用地分布图，山洪淹没区域的各水文站历史观测资料，山洪淹没区域的水系选定河道各断面水位流量关系、水库和湖泊水位库容关系，山洪淹没区域选定的河网水系分布，河道水下地形图，河道边界上的水位和流量资料；
(2)采用无规则非结构网格针对复杂计算区域和不规则地形进行网格剖分，使计算区域离散化；
(3)逐时段用有限体积法对每一单元网格进行水量、动量平衡计算，用黎曼近似解计算跨单元的水量和动量的法向数值通量，计算出连续时段不同单元网格的水量Q和坐标轴X方向的动量U和坐标轴Y方向的动量V；
(4)根据步骤(3)得出的结果，进一步计算出山洪抵达某一单元网格所需的时段，获得每个单元网格的洪水到达时间，依据对单元网格水量的统计值和所经历的时段获得某一单元网格的洪水淹没时间，同理根据步骤(3)的计算结果计算得出最大洪水淹没水深、最大洪水水位和最大洪水流速；
(5)计算每个单元网格洪水到达时间、洪水淹没时间、最大洪水淹没水深，最大的洪水水位和最大洪水流速，分别生成等值线，得到洪水到达时间等值图、洪水淹没时间等值图、最大洪水淹没水深等值图、最大的洪水水位等值图和最大洪水流速等值图。
优选地，网格剖分采用四边形网格、三角形网格或四边形网格与三角形网格的混合，包括以下步骤：
21)对山洪淹没区域内的水工建筑物和地貌数据进行概化处理，如滩地、荒地、农田、林地、道路、堤防、民区、河渠、湖泊、水库等各种复杂的建筑和地貌，形成计算网格的边界控制点，进而确定模拟区域的内部边界，内部边界由一个或多个关键控制点构成；
22)计算网格的内部控制点，包括模拟区域水面之外地形的散点和水下地形的散点，水面之外的地形控制点为对地形数据等间距提取出设定精度的散点，形成计算网格的内部控制点；水下地形利用河道大断面测绘数据提取设定精度的散点，沿河道行进方向控制点间距与大断面测绘间距一致，形成计算网格的内部控制点，用来计算山洪水流河道内的演进过程；
23)定义上游、下游、左岸、右岸控制边界点，边界点即为三角形网格的边界计算控制点；
24)在模拟区域中根据控制点剖分三角形网格，；
25)通过边界尺度控制网格大小及网格密度，即对洪水演进模拟重点区域加密网格密度控制边界点，在设定边界控制点的同时标识各类边界的类型，包括水位边界、流量边界、水位流量关系边界和岸边界，对于水工结构边界，根据类型特点单独设置边界，包括给定闸门泄流能力公式和溃口发展过程。
优选地，加密网格密度控制边界点的方式包括：选择边界进行均匀等距剖分或渐进加密剖分，或者以上方式的组合形式。
优选地，对山洪淹没区域内地形复杂的区域利用非结构三角形网格进行加密剖分，实现特定区域的内部精细网格剖分，具体步骤为：
41)河道以外区域通过数字高程模型DEM等高线，提取出设定精度的散点；
42)获取河道内的水下地形散点，根据水文站河道断面的水位和流量信息，即断面关系曲线和实测数据，获得水深，通过水位-水深＝水下地形，计算出断面某一点的水下地形高程，逐断面计算获得河道每个断面上多组水下地形散点数据，对于缺失的水文站断面，依次进行推算；
43)将河道外和河道内的地形散点合并，构成三角形网格构建所需的地形散点，计算三角形网格并通过插值生成高精度三角形网格，即采用邻点相连的方法，将每个中心点链接成三角形网格，对于边界上的点，将边界点和邻近的中心点，组成三角形网格，边界上的点的值根据邻近的中心点的值取平均值；
44)三角形网格生成后，通过插值获取三角形网格单元中心的高程。
优选地，对每一单元网格进行水量、动量平衡计算的步骤为：
51)山洪动态演进仿真模型的基本控制方程为：
Figure BDA0001556991010000041
式中q＝[h,hu,hv]T为守恒物理向量,f(q)＝[hu,hu2+gh2/2,huv]T为坐标轴x向的通量向量，g(q)＝[hv,huv,hv2+gh2/2]T为坐标轴y向的通量向量；h为水深，u和v分别为坐标轴x和y向垂线平均流速分量，g为重力加速度；源汇项b(q)为
b(q)＝[qw,gh(s0x-sfx)+qwu,gh(s0y-sfy)+qwv]
其中s0x和sfx分别为坐标轴x向的河底坡度及摩阻坡度，s0y和sfy分别为坐标轴y向的河底坡度及摩阻坡度，qw为单位时间内的净雨深；所述仿真模型中摩阻坡度由曼宁公式估算；
52)应用散度定理对所述山洪动态演进仿真模型基本控制方程在任意单元Ω上进行积分离散，得到有限体积法的基本方程
Figure BDA0001556991010000054
其中qt为方程组的时变项，通过守恒物理向量q对t求导求得，n为单元边
Figure BDA0001556991010000055
外法向单位向量，dω和dL分别为面积分微元和线积分微元，F(q)·n为法向数值通量Fn(q)，F(q)＝[f(q),g(q)]T；向量q为单元平均值，对于一阶精度离散假定为常数；上述方程离散化求得到：
Figure BDA0001556991010000051
其中b*(q)＝(A·b1,A·b2,A·b3)T，式中A为单元Ω的面积，b1＝qw，b2＝gh(s0x-sfx)+qwu，b3＝gh(s0y-sfy)+qwv，m为单元边总数，Lj为单元中第j边的长度；b*(q)为源汇项，单元边法向数值通量
Figure BDA0001556991010000052
简记为Fn(q)，定义如下：
Fn(q)＝cosΦ·f(q)+sinΦ·g(q)
根据f(q)和g(q)的坐标变换旋转不变性，即满足
Figure BDA0001556991010000053
得到
Figure BDA0001556991010000061
其中，Φ为单元面法向向量n与坐标轴x轴的夹角，夹角由x轴起逆时针计量；T(Φ)和T(Φ)-1分别为坐标旋转变换矩阵及其逆阵，表达式为：
Figure BDA0001556991010000062
可得
Figure BDA0001556991010000063
式中
Figure BDA0001556991010000064
向量
Figure BDA0001556991010000065
是向量q在法向上的投影，相应的流速分量分别为法向和切向；
53)采用黎曼问题计算法向通量
Figure BDA0001556991010000066
求出下一时刻的向量q。
优选地，利用下式计算跨单元的水量、动量的法向数值通量：
Figure BDA0001556991010000067
满足：
Figure BDA0001556991010000068
Figure BDA0001556991010000069
为单元界面上法向量n的坐标，原点位于单元边中点，轴向与外法向一致，
Figure BDA00015569910100000610
为该局部坐标原点处的外法向通量，
Figure BDA00015569910100000611
通过守恒物理向量
Figure BDA00015569910100000620
对时间t求导求得，
Figure BDA00015569910100000612
为
Figure BDA00015569910100000613
对
Figure BDA00015569910100000614
求导得到；
Figure BDA00015569910100000615
和
Figure BDA00015569910100000616
分别为向量
Figure BDA00015569910100000617
在单元界面左右的状态，模型约定计算单元为左边，相邻单元则为右边；t＝0时的初始状态已知，通过解算黎曼问题，得到原点位于
Figure BDA00015569910100000618
时间为t＝0+的外法向数值通量，记为fLR(qL,qR)，进而得到法向通量
Figure BDA00015569910100000619
优选地，法向通量
Figure BDA0001556991010000071
的计算方法包括：由公共边两侧单元通量取平均值计算通量，
Figure BDA0001556991010000072
或由两侧单元物理守恒量的均值计算通量，
Figure BDA0001556991010000073
利用单调性格式计算通量，包括全变差缩小格式TVD和通量输运校正格式FCT；利用基于特征理论并具有逆风性的黎曼近似解计算通量，包括通量向量分裂格式FVS、通量差分裂格式FDS和Osher格式。
与现有技术相比，本发明所达到的有益效果是：本发明针对复杂计算域和不规则地形的山洪易发山丘区，提供了一种基于无规则非结构网格模型的山洪动态演进仿真方法，可实现洪水淹没过程的动态模拟和仿真，实时模拟和计算山洪淹没范围、淹没水深、流场等洪水要素，实现洪水淹没过程的可视化仿真。
附图说明
图1为本发明方法中采用的山洪仿真模型计算流程图；
图2为本发明的山洪仿真模拟区域选取示意图；
图3为本发明的无规则非结构网格剖分和插值图；
图4为本发明的时段山洪水深等值线图。
具体实施方式
下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
图1为本发明方法中采用的山洪仿真模型计算流程图，本发明的一种基于无规则非结构网格模型的山洪动态仿真方法，具体包括以下步骤：
(1)选定模拟区域，山洪易发的典型区域，或易发生垮堤的土堤河段，输入资料数据包括：山洪淹没区域的气候、水文、地形、地貌、地质及水工建筑物基本概况资料，山洪淹没区域的水上数字高程模型DEM地形数据和河道地形数据，山洪淹没区域的分类用地分布图，山洪淹没区域的各水文站历史观测资料，山洪淹没区域的水系选定河道各断面水位流量关系、水库和湖泊水位库容关系，山洪淹没区域选定的河网水系分布，河道水下地形图，河道边界上的水位和流量资料；
(2)采用无规则非结构网格针对复杂计算区域和不规则地形进行网格剖分，使计算区域离散化；
(3)逐时段用有限体积法对每一单元网格进行水量、动量平衡计算，用黎曼近似解计算跨单元的水量和动量的法向数值通量，计算出连续时段不同单元网格的水量Q和坐标轴X方向的动量U和坐标轴Y方向的动量V；
(4)根据步骤(3)得出的结果，进一步计算出山洪抵达某一单元网格所需的时段，获得每个单元网格的洪水到达时间，依据对单元网格水量的统计值和所经历的时段获得某一单元网格的洪水淹没时间，同理根据步骤(3)的计算结果计算得出最大洪水淹没水深、最大洪水水位和最大洪水流速；
(5)计算每个单元网格洪水到达时间、洪水淹没时间、最大洪水淹没水深，最大的洪水水位和最大洪水流速，分别生成等值线，得到洪水到达时间等值图、洪水淹没时间等值图、最大洪水淹没水深等值图、最大的洪水水位等值图和最大洪水流速等值图，图4为本发明的时段山洪水深等值线图。
优选地，图3为本发明的无规则非结构网格剖分和插值图，网格剖分采用四边形网格、三角形网格或四边形网格与三角形网格的混合，包括以下步骤：
21)对山洪淹没区域内的水工建筑物和地貌数据进行概化处理，如滩地、荒地、农田、林地、道路、堤防、民区、河渠、湖泊、水库等各种复杂的建筑和地貌，形成计算网格的边界控制点，进而确定模拟区域的内部边界，内部边界由一个或多个关键控制点构成，如图2为本发明的山洪仿真模拟区域选取示意图；
22)计算网格的内部控制点，包括模拟区域水面之外地形的散点和水下地形的散点，水面之外的地形控制点为对高精度地形数据如1:500或者1:5000等间距提取出设定精度的散点，形成计算网格的内部控制点；水下地形利用河道大断面测绘数据提取设定精度的散点，沿河道行进方向控制点间距与大断面测绘间距一致，形成计算网格的内部控制点，用来计算山洪水流河道内的演进过程；
23)定义上游、下游、左岸、右岸控制边界点，边界点即为三角形网格的边界计算控制点；
24)在模拟区域中根据控制点剖分三角形网格，；
25)通过边界尺度控制网格大小及网格密度，即对洪水演进模拟重点区域加密网格密度控制边界点，在设定边界控制点的同时标识各类边界的类型，包括水位边界、流量边界、水位流量关系边界和岸边界，对于水工结构边界，根据类型特点单独设置边界，包括给定闸门泄流能力公式和溃口发展过程。
优选地，加密网格密度控制边界点的方式包括：选择边界进行均匀等距剖分或渐进加密剖分，或者以上方式的组合形式。
优选地，对山洪淹没区域内地形复杂的区域利用非结构三角形网格进行加密剖分，实现特定区域的内部精细网格剖分，具体步骤为：
41)河道以外区域通过数字高程模型DEM等高线，提取出设定精度的散点；
42)获取河道内的水下地形散点，根据水文站河道断面的水位和流量信息，即断面关系曲线和实测数据，获得水深，通过水位-水深＝水下地形，计算出断面某一点的水下地形高程，逐断面计算获得河道每个断面上多组水下地形散点数据，对于缺失的水文站断面，依次进行推算；
43)将河道外和河道内的地形散点合并，构成三角形网格构建所需的地形散点，计算三角形网格并通过插值生成高精度三角形网格，即采用邻点相连的方法，将每个中心点链接成三角形网格，对于边界上的点，将边界点和邻近的中心点，组成三角形网格，边界上的点的值根据邻近的中心点的值取平均值；
44)三角形网格生成后，通过插值获取三角形网格单元中心的高程。
优选地，对每一单元网格进行水量、动量平衡计算的步骤为：
51)山洪动态演进仿真模型的基本控制方程为：
Figure BDA0001556991010000101
式中q＝[h,hu,hv]T为守恒物理向量,f(q)＝[hu,hu2+gh2/2,huv]T为坐标轴x向的通量向量，g(q)＝[hv,huv,hv2+gh2/2]T为坐标轴y向的通量向量；h为水深，u和v分别为坐标轴x和y向垂线平均流速分量，g为重力加速度；源汇项b(q)为
b(q)＝[qw,gh(s0x-sfx)+qwu,gh(s0y-sfy)+qwv]
其中s0x和sfx分别为坐标轴x向的河底坡度及摩阻坡度，s0y和sfy分别为坐标轴y向的河底坡度及摩阻坡度，qw为单位时间内的净雨深；所述仿真模型中摩阻坡度由曼宁公式估算；
52)应用散度定理对所述山洪动态演进仿真模型基本控制方程在任意单元Ω上进行积分离散，得到有限体积法的基本方程
Figure BDA0001556991010000102
其中qt为方程组的时变项，通过守恒物理向量q对t求导求得，n为单元边
Figure BDA0001556991010000111
外法向单位向量，dω和dL分别为面积分微元和线积分微元，F(q)·n为法向数值通量Fn(q)，F(q)＝[f(q),g(q)]T；向量q为单元平均值，对于一阶精度离散假定为常数；上述方程离散化求得到：
Figure BDA0001556991010000112
其中b*(q)＝(A·b1,A·b2,A·b3)T，式中A为单元Ω的面积，b1＝qw，b2＝gh(s0x-sfx)+qwu，b3＝gh(s0y-sfy)+qwv，m为单元边总数，Lj为单元中第j边的长度；b*(q)为源汇项，单元边法向数值通量
Figure BDA0001556991010000119
(q)简记为Fn(q)，定义如下：
Fn(q)＝cosΦ·f(q)+sinΦ·g(q)
根据f(q)和g(q)的坐标变换旋转不变性，即满足
Figure BDA0001556991010000113
得到
Figure BDA0001556991010000114
其中，Φ为单元面法向向量n与坐标轴x轴的夹角，夹角由x轴起逆时针计量；T(Φ)和T(Φ)-1分别为坐标旋转变换矩阵及其逆阵，表达式为：
Figure BDA0001556991010000115
可得
Figure BDA0001556991010000116
式中
Figure BDA0001556991010000117
向量
Figure BDA0001556991010000118
是向量q在法向上的投影，相应的流速分量分别为法向和切向；
53)采用黎曼问题计算法向通量
Figure BDA0001556991010000121
求出下一时刻的向量q。
优选地，利用下式计算跨单元的水量、动量的法向数值通量：
Figure BDA0001556991010000122
满足：
Figure BDA0001556991010000123
Figure BDA0001556991010000124
为单元界面上法向量n的坐标，原点位于单元边中点，轴向与外法向一致，
Figure BDA0001556991010000125
为该局部坐标原点处的外法向通量，
Figure BDA0001556991010000126
通过守恒物理向量
Figure BDA0001556991010000127
对时间t求导求得，
Figure BDA0001556991010000128
为
Figure BDA0001556991010000129
对
Figure BDA00015569910100001210
求导得到；
Figure BDA00015569910100001211
和
Figure BDA00015569910100001212
分别为向量
Figure BDA00015569910100001213
在单元界面左右的状态，模型约定计算单元为左边，相邻单元则为右边；t＝0时的初始状态已知，通过解算黎曼问题，得到原点位于
Figure BDA00015569910100001214
时间为t＝0+的外法向数值通量，记为fLR(qL,qR)，进而得到法向通量
Figure BDA00015569910100001215
优选地，法向通量
Figure BDA00015569910100001216
的计算方法包括：由公共边两侧单元通量取平均值计算通量，
Figure BDA00015569910100001217
或由两侧单元物理守恒量的均值计算通量，
Figure BDA00015569910100001218
利用单调性格式计算通量，包括全变差缩小格式TVD和通量输运校正格式FCT；利用基于特征理论并具有逆风性的黎曼近似解计算通量，包括通量向量分裂格式FVS、通量差分裂格式FDS和Osher格式。

Claims (5)

1.一种基于无规则非结构网格模型的山洪动态演进仿真方法，其特征在于，包括以下步骤：
(1)选定模拟区域，输入资料数据包括：山洪淹没区域的气候、水文、地形、地貌、地质及水工建筑物基本概况资料，山洪淹没区域的水上数字高程模型DEM地形数据和河道地形数据，山洪淹没区域的分类用地分布图，山洪淹没区域的各水文站历史观测资料，山洪淹没区域的水系选定河道各断面水位流量关系、水库和湖泊水位库容关系，山洪淹没区域选定的河网水系分布，河道水下地形图，河道边界上的水位和流量资料；
(2)采用无规则非结构网格针对复杂计算区域和不规则地形进行网格剖分，使计算区域离散化，其中，网格剖分采用四边形网格、三角形网格或四边形网格与三角形网格的混合，包括以下步骤：
21)对山洪淹没区域内的水工建筑物和地貌数据进行概化处理，形成计算网格的边界控制点，进而确定模拟区域的内部边界；
22)计算网格的内部控制点，包括模拟区域水面之外地形的散点和水下地形的散点，水面之外的地形控制点为对地形数据等间距提取散点，形成计算网格的内部控制点；水下地形利用河道大断面测绘数据提取设定精度的散点，沿河道行进方向控制点间距与大断面测绘间距一致，形成计算网格的内部控制点；
23)定义上游、下游、左岸、右岸控制边界点，边界点即为三角形网格的边界计算控制点；
24)在模拟区域中根据控制点剖分三角形网格；
25)通过边界尺度控制网格大小及网格密度，即对洪水演进模拟重点区域加密网格密度控制边界点，在设定边界控制点的同时标识各类边界的类型，包括水位边界、流量边界、水位流量关系边界和岸边界，对于水工结构边界，根据类型特点单独设置边界；其中，对山洪淹没区域内地形复杂的区域利用非结构三角形网格进行加密剖分，实现特定区域的内部精细网格剖分，具体步骤为：
41)河道以外区域通过数字高程模型DEM等高线，提取出设定精度的散点；
42)获取河道内的水下地形散点，根据水文站河道断面的水位和流量信息，即断面关系曲线和实测数据，获得水深，通过水位-水深＝水下地形，计算出断面某一点的水下地形高程，逐断面计算获得河道每个断面上多组水下地形散点数据，对于缺失的水文站断面，依次进行推算；
43)将河道外和河道内的地形散点合并，构成三角形网格构建所需的地形散点，计算三角形网格并通过插值生成高精度三角形网格，即采用邻点相连的方法，将每个中心点链接成三角形网格，对于边界上的点，将边界点和邻近的中心点，组成三角形网格，边界上的点的值根据邻近的中心点的值取平均值；
44)三角形网格生成后，通过插值获取三角形网格单元中心的高程；
(3)逐时段用有限体积法对每一单元网格进行水量、动量平衡计算，用黎曼近似解计算跨单元的水量和动量的法向数值通量，计算出连续时段不同单元网格的水量Q和坐标轴X方向的动量U和坐标轴Y方向的动量V；
(4)根据步骤(3)得出的结果，进一步计算出山洪抵达某一单元网格所需的时段，获得每个单元网格的洪水到达时间，依据对单元网格水量的统计值和所经历的时段获得某一单元网格的洪水淹没时间，同理根据步骤(3)的计算结果计算得出最大洪水淹没水深、最大洪水水位和最大洪水流速；
(5)计算每个单元网格洪水到达时间、洪水淹没时间、最大洪水淹没水深，最大的洪水水位和最大洪水流速，分别生成等值线，得到洪水到达时间等值图、洪水淹没时间等值图、最大洪水淹没水深等值图、最大的洪水水位等值图和最大洪水流速等值图。
2.根据权利要求1所述的一种基于无规则非结构网格模型的山洪动态演进仿真方法，其特征在于，加密网格密度控制边界点的方式包括：选择边界进行均匀等距剖分或渐进加密剖分，或者以上方式的组合形式。
3.根据权利要求1所述的一种基于无规则非结构网格模型的山洪动态演进仿真方法，其特征在于，在所述步骤(3)中，对每一单元网格进行水量、动量平衡计算的步骤为：
51)山洪动态演进仿真模型的基本控制方程为：
Figure FDA0003073922860000031
式中q＝[h,hu,hv]T为守恒物理向量,f(q)＝[hu,hu2+gh2/2,huv]T为坐标轴x向的通量向量，g(q)＝[hv,huv,hv2+gh2/2]T为坐标轴y向的通量向量；h为水深，u和v分别为坐标轴x和y向垂线平均流速分量，g为重力加速度；源汇项b(q)为
b(q)＝[qw,gh(s0x-sfx)+qwu,gh(s0y-sfy)+qwv]
其中s0x和sfx分别为坐标轴x向的河底坡度及摩阻坡度，s0y和sfy分别为坐标轴y向的河底坡度及摩阻坡度，qw为单位时间内的净雨深；所述仿真模型中摩阻坡度由曼宁公式估算；
52)应用散度定理对所述山洪动态演进仿真模型基本控制方程在任意单元Ω上进行积分离散，得到有限体积法的基本方程
Figure FDA0003073922860000032
其中qt为方程组的时变项，通过守恒物理向量q对t求导求得，n为单元边
Figure FDA0003073922860000041
外法向单位向量，dω和dL分别为面积分微元和线积分微元，F(q)·n为法向数值通量Fn(q)，F(q)＝[f(q),g(q)]T；守恒物理向量q为单元平均值，对于一阶精度离散假定为常数；上述方程离散化求得到：
Figure FDA0003073922860000042
其中，b(q)＝(A·b1,A·b2,A·b3)T,式中A为单元Ω的面积，b1＝qw，b2＝gh(s0x-sfx)+qwu，b3＝gh(s0y-sfy)+qwv，m为单元边总数，Lj为单元中第j边的长度；b(q)为源汇项，
Figure FDA0003073922860000043
简记为Fn(q)，定义如下：
Fn(q)＝cosΦ·f(q)+sinΦ·g(q)
根据f(q)和g(q)的坐标变换旋转不变性，即满足
Figure FDA0003073922860000044
得到
Figure FDA0003073922860000045
其中，Φ为单元边
Figure FDA0003073922860000046
外法向单位向量n与坐标轴x轴的夹角，夹角由x轴起逆时针计量；T(Φ)和T(Φ)-1分别为坐标旋转变换矩阵及其逆阵，表达式为：
Figure FDA0003073922860000047
可得
Figure FDA0003073922860000048
式中
Figure FDA0003073922860000049
向量
Figure FDA00030739228600000410
是守恒物理向量q在法向上的投影，相应的流速分量分别为法向和切向；
53)采用黎曼问题计算向量
Figure FDA0003073922860000051
的法向通量
Figure FDA0003073922860000052
求出下一时刻的守恒物理向量q。
4.根据权利要求3所述的一种基于无规则非结构网格模型的山洪动态演进仿真方法，其特征在于，在所述步骤(3)中，利用下式计算跨单元的水量、动量的法向数值通量：
Figure FDA0003073922860000053
满足：
Figure FDA0003073922860000054
Figure FDA0003073922860000055
为单元界面上法向量n的坐标，原点位于单元边中点，轴向与外法向一致，
Figure FDA0003073922860000056
通过向量
Figure FDA0003073922860000057
对时间t求导求得，
Figure FDA0003073922860000058
由
Figure FDA0003073922860000059
对
Figure FDA00030739228600000510
求导得到；
Figure FDA00030739228600000511
和
Figure FDA00030739228600000512
分别为向量
Figure FDA00030739228600000513
在单元界面左右的状态，模型约定计算单元为左边，相邻单元则为右边；t＝0时的初始状态已知，通过解算黎曼问题，得到原点位于
Figure FDA00030739228600000514
时间为t＝0的外法向数值通量，记为fLR(qL,qR)，进而得到向量
Figure FDA00030739228600000515
的法向通量
Figure FDA00030739228600000516
5.根据权利要求4所述的一种基于无规则非结构网格模型的山洪动态演进仿真方法，其特征在于，法向通量
Figure FDA00030739228600000517
的计算方法包括：由公共边两侧单元通量取平均值计算通量，
Figure FDA00030739228600000518
或由两侧单元物理守恒量的均值计算通量，
Figure FDA00030739228600000519
利用单调性格式计算通量，包括全变差缩小格式TVD和通量输运校正格式FCT；利用基于特征理论并具有逆风性的黎曼近似解计算通量，包括通量向量分裂格式FVS、通量差分裂格式FDS和Osher格式。
CN201810067029.XA 2018-01-24 2018-01-24 一种基于无规则非结构网格模型的山洪动态仿真方法 Active CN108021780B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810067029.XA CN108021780B (zh) 2018-01-24 2018-01-24 一种基于无规则非结构网格模型的山洪动态仿真方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810067029.XA CN108021780B (zh) 2018-01-24 2018-01-24 一种基于无规则非结构网格模型的山洪动态仿真方法

Publications (2)

Publication Number Publication Date
CN108021780A CN108021780A (zh) 2018-05-11
CN108021780B true CN108021780B (zh) 2021-08-13

Family

ID=62074702

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810067029.XA Active CN108021780B (zh) 2018-01-24 2018-01-24 一种基于无规则非结构网格模型的山洪动态仿真方法

Country Status (1)

Country Link
CN (1) CN108021780B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109543275B (zh) * 2018-11-15 2019-07-09 中国水利水电科学研究院 一种城区地表径流二维数值模拟方法
CN109523078B (zh) * 2018-11-16 2020-08-18 黄河勘测规划设计研究院有限公司 一种洪水风险图的优化方法及系统
CN111191372B (zh) * 2020-01-01 2021-12-21 浙江大学 一种基于二维水动力学的大规模洪水场景模拟预警交互方法和系统
CN111241757A (zh) * 2020-01-10 2020-06-05 中核第四研究设计工程有限公司 基于计算流体力学的铀尾矿库溃坝三维数值模拟方法
CN112257313B (zh) * 2020-10-21 2024-05-14 西安理工大学 一种基于gpu加速的污染物输移高分辨率数值模拟方法
CN114170397B (zh) * 2022-02-09 2022-04-29 四川省安全科学技术研究院 基于真实地形的不规则离散元仿真模型快速数学建模方法
CN116150548B (zh) * 2023-04-17 2023-07-21 云南省水利水电科学研究院 一种河道洪水淹没范围计算方法
CN116522818B (zh) * 2023-05-09 2023-12-19 中国水利水电科学研究院 一种大坡度地形条件下的干旱区潜水位模拟方法
CN117173369B (zh) * 2023-11-02 2023-12-29 华中科技大学 基于WebGL的三维洪水演进模拟方法及系统

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102496168A (zh) * 2011-11-22 2012-06-13 南京大学 一种用于河道水文数值模拟的复杂河道网格化方法
KR20140125136A (ko) * 2013-04-18 2014-10-28 창원대학교 산학협력단 격자생성 프로그램을 이용한 지형 분석방법
CN104851360A (zh) * 2014-02-14 2015-08-19 杭州贵仁科技有限公司 一种洪水风险图的生成方法和系统
CN105844709A (zh) * 2016-03-25 2016-08-10 中国水利水电科学研究院 复杂河道地形流域洪水演进虚拟仿真的淹没线追踪方法
CN106599140A (zh) * 2016-12-06 2017-04-26 河海大学 一种基于gis的洪水风险要素的快速处理方法
CN107239657A (zh) * 2017-05-31 2017-10-10 中国水利水电科学研究院 一种面向对象的水动力学建模要素管理方法
CN108062453A (zh) * 2018-01-12 2018-05-22 河南省水利勘测设计研究有限公司 水利信息化系统洪水高效模拟与高逼真可视化动态展示方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8768663B2 (en) * 2010-10-08 2014-07-01 The United States Of America, As Represented By The Secretary Of The Navy Automated method and system for predicting high resolution tidal heights and currents in coastal zones
US8676555B2 (en) * 2010-10-26 2014-03-18 The United States Of America, As Represented By The Secretary Of The Navy Tool for rapid configuration of a river model using imagery-based information

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102496168A (zh) * 2011-11-22 2012-06-13 南京大学 一种用于河道水文数值模拟的复杂河道网格化方法
KR20140125136A (ko) * 2013-04-18 2014-10-28 창원대학교 산학협력단 격자생성 프로그램을 이용한 지형 분석방법
CN104851360A (zh) * 2014-02-14 2015-08-19 杭州贵仁科技有限公司 一种洪水风险图的生成方法和系统
CN105844709A (zh) * 2016-03-25 2016-08-10 中国水利水电科学研究院 复杂河道地形流域洪水演进虚拟仿真的淹没线追踪方法
CN106599140A (zh) * 2016-12-06 2017-04-26 河海大学 一种基于gis的洪水风险要素的快速处理方法
CN107239657A (zh) * 2017-05-31 2017-10-10 中国水利水电科学研究院 一种面向对象的水动力学建模要素管理方法
CN108062453A (zh) * 2018-01-12 2018-05-22 河南省水利勘测设计研究有限公司 水利信息化系统洪水高效模拟与高逼真可视化动态展示方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
MAST-2D diffusive model for flood prediction on domains with triangular Delaunay unstructured meshes;C. Aricò et al.;《Advances in Water Resources》;20110905(第34期);全文 *
基于GIS的二维溃坝洪水演进数值模拟及可视化系统研制;王晓航;《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》;20071115(第5期);第17-35，41-53页 *
王晓航.基于GIS的二维溃坝洪水演进数值模拟及可视化系统研制.《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》.2007,(第5期), *

Also Published As

Publication number Publication date
CN108021780A (zh) 2018-05-11

Similar Documents

Publication Publication Date Title
CN108021780B (zh) 一种基于无规则非结构网格模型的山洪动态仿真方法
CN107506566B (zh) 一种泥石流动力学数值模拟分析方法及系统
Fettweis et al. The mud deposits and the high turbidity in the Belgian–Dutch coastal zone, southern bight of the North Sea
Worni et al. Coupling glacial lake impact, dam breach, and flood processes: A modeling perspective
Mizukami et al. mizuRoute version 1: a river network routing tool for a continental domain water resources applications
CN107451372B (zh) 一种运动波与动力波相结合的山区洪水过程数值模拟方法
Zhang et al. Numerical simulation and analysis of saltwater intrusion lengths in the Pearl River Delta, China
CN112784502B (zh) 一种水文水力学动态双向耦合的洪水预测方法
Wang et al. Long-term evolution in the location, propagation, and magnitude of the tidal shear front off the Yellow River Mouth
Poretti et al. An approach for flood hazard modelling and mapping in the medium Valtellina
Jeannot et al. A low-dimensional integrated subsurface hydrological model coupled with 2-D overland flow: Application to a restored fluvial hydrosystem (Upper Rhine River–France)
Comer et al. Development of high-resolution multi-scale modelling system for simulation of coastal-fluvial urban flooding
Roushangar et al. Linear and non-linear approaches to predict the Darcy-Weisbach friction factor of overland flow using the extreme learning machine approach
Zhu et al. Morphologic modelling of tidal inlet on a barrier-lagoon coast: Case study of the Laolonggou tidal inlet in the Bohai Bay
Busaman et al. Dynamically adaptive tree grid modeling for simulation and visualization of rainwater overland flow
Krylenko et al. Analysis of the impact of hydrotechnical construction on the Amur river near Blagoveshchensk and Heihe cities using a two-dimensional hydrodynamic model
Gergeľová et al. Hydrodynamic modeling and GIS tools applied in urban areas.
Thakur et al. Exploring CCHE2D and its sediment modelling capabilities
Martin et al. MOD_FreeSurf2D: A MATLAB surface fluid flow model for rivers and streams
Wang Numerical improvements for large-scale flood simulation
David et al. Basics for hydraulic modelling of flood runoff using advanced hydroinformatic tools
Ding et al. Integrated coastal process modeling and impact assessment of flooding and sedimentation in coasts and estuaries
Isaac et al. Experimental and numerical simulations of the hydraulic flushing of reservoir sedimentation for planning and design of a hydropower project
Solaimani Flood forecasting based on GIS and hydraulic model
Jetten Catchment-scale multi-process modeling with local time stepping

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant