CN110765694B - 一种基于简化型浅水方程组的城市地表水流数值模拟方法 - Google Patents
一种基于简化型浅水方程组的城市地表水流数值模拟方法 Download PDFInfo
- Publication number
- CN110765694B CN110765694B CN201911147817.0A CN201911147817A CN110765694B CN 110765694 B CN110765694 B CN 110765694B CN 201911147817 A CN201911147817 A CN 201911147817A CN 110765694 B CN110765694 B CN 110765694B
- Authority
- CN
- China
- Prior art keywords
- water
- elevation
- unit
- interface
- equation set
- 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
Links
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 title claims abstract description 138
- 238000004088 simulation Methods 0.000 title claims abstract description 54
- 238000000034 method Methods 0.000 title claims abstract description 51
- 239000002352 surface water Substances 0.000 title claims abstract description 39
- 230000004907 flux Effects 0.000 claims abstract description 38
- 238000004364 calculation method Methods 0.000 claims abstract description 29
- 238000011160 research Methods 0.000 claims abstract description 23
- 230000008569 process Effects 0.000 claims abstract description 13
- 239000006185 dispersion Substances 0.000 claims abstract description 12
- 210000004027 cell Anatomy 0.000 claims description 55
- 239000011159 matrix material Substances 0.000 claims description 32
- 239000013598 vector Substances 0.000 claims description 15
- 230000001133 acceleration Effects 0.000 claims description 11
- 238000012545 processing Methods 0.000 claims description 8
- 230000002706 hydrostatic effect Effects 0.000 claims description 6
- 238000012937 correction Methods 0.000 claims description 5
- 230000000694 effects Effects 0.000 claims description 4
- 230000005484 gravity Effects 0.000 claims description 4
- 238000012986 modification Methods 0.000 claims description 3
- 230000004048 modification Effects 0.000 claims description 3
- 210000004457 myocytus nodalis Anatomy 0.000 claims description 2
- 230000033001 locomotion Effects 0.000 abstract description 5
- 230000003068 static effect Effects 0.000 description 4
- 230000008859 change Effects 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 239000012530 fluid Substances 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- YTAHJIFKAKIKAV-XNMGPUDCSA-N [(1R)-3-morpholin-4-yl-1-phenylpropyl] N-[(3S)-2-oxo-5-phenyl-1,3-dihydro-1,4-benzodiazepin-3-yl]carbamate Chemical compound O=C1[C@H](N=C(C2=C(N1)C=CC=C2)C1=CC=CC=C1)NC(O[C@H](CCN1CCOCC1)C1=CC=CC=C1)=O YTAHJIFKAKIKAV-XNMGPUDCSA-N 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000005293 physical law Methods 0.000 description 1
- 230000008092 positive effect Effects 0.000 description 1
- 238000000746 purification Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供了一种基于简化型浅水方程组的城市地表水流数值模拟方法。本发明首先获取城市地区的地形、地貌数据后,对研究区域用三角形非结构网格进行离散,考虑城市区域地表水流的运动特征,采用忽略对流项的简化型浅水方程组,并将其Godunov型的离散格式作为数值模拟的控制方程;采用Roe格式黎曼求解器对单元界面处形成的黎曼问题进行求解,并通过一种新的水面重构方法和隐式方法分别对底坡源项以及摩擦源项进行空间离散,构建了一种基于简化型浅水方程组的城市地表水流数值模型。本发明通过对浅水方程组进行简化,大量减少求解过程中复杂的通量计算,提高了城市地表水流数值模拟的效率,为城市地表水流数值模拟提供了一种新的途径。
Description
技术领域
本发明涉及计算流体力学领域,具体涉及一种基于简化型浅水方程组的城市地表水流数值模拟方法。
背景技术
近年来,极端降雨事件频发,城市暴内涝现象呈现快速增长趋势,洪涝灾害带来的影响已成为制约当地社会经济发展的重要因素。面对日益严峻的洪涝灾害,如何完善城市洪涝预警预报系统,减轻由洪涝灾害造成的社会经济损失,已经成为目前研究的热点和难点,而数值模型提供了一种以较低成本模拟城市洪水演进过程的有效解决方案。
城市空间内分布着密集的不透水建筑物,同时城区建筑物轮廓往往是不规则的,这对城市地表水流数值模拟的计算效率带来严峻挑战。相对于自然流域,城市区域由于存在着密集的不透水区域,城市雨洪模型的网格数量将变得极为庞大,同时在数值模拟过程中,基于完整二维浅水方程组的数值模型虽然能更好地体现洪水运动的物理特征,但通常需要复杂的通量计算来提供稳定精确的数值解(Hongbin Zhang,2014:Non-negativedepth reconstruction for a two-dimensionalpartial inertial inundation model),这成为限制城市洪涝预警预报系统发展的重要因素。
发明内容
为了提高城市地表水流数值模拟的效率,本发明摒弃了现有城市雨洪模型中常用的完整浅水方程组,提出了一种简化型浅水方程组的有限体积解法,用于城市地表水流的数值模拟,通过减少每个网格单元上的通量计算以提高城市地表水流数值模拟的效率,应用于城市洪涝预警预报系统以提高其效率。
本发明的目的至少通过以下技术方案之一实现。
一种基于简化型浅水方程组的城市地表水流数值模拟方法,包括以下步骤:
步骤一、获取研究区域的地理数字高程信息、土地利用类型遥感影像信息以及不透水建筑物的几何轮廓信息;
步骤二、采用三角形非结构网格对研究区域进行空间离散,其中不透水建筑物部分不进行网格剖分,将步骤一获得的研究区域的数字高程信息及土地利用类型遥感影像信息在网格上插值得到每个网格的高程值及糙率值;
步骤三、对完整的二维浅水方程组进行简化,推导简化型浅水方程组的守恒格式及非线性矩阵形式,得到用于构建数值模型的控制方程,并基于所述控制方程对城市地表水流进行数值模拟;
步骤四、为了适用城市地表特征,采用有限体积法对守恒格式简化型浅水方程组的控制方程组在三角形非结构网格上进行离散,同时确定各物理变量,包括水深、流速、动量和质量在网格中的位置,最终得到简化型浅水方程组的离散后的积分格式;
步骤五、针对简化型浅水方程数值模型单元间的流量通量计算中存在的黎曼问题,采用Roe格式黎曼求解器,构造相应的近似雅克比矩阵以及由其特征值和特征向量组成的相关系数矩阵,获得单元界面处各物理变量的黎曼状态,以解决单元界面处左右物理量不一致的问题;
步骤六、通过一种新的水面重构方法,修正黎曼状态量的表达式,并对简化型浅水方程组的底坡源项进行离散,得到底坡源项的迎风离散格式;
步骤七、针对城市地表水流数值模拟中的干湿水深动边界问题进行处理,得到修正后的底坡源项的迎风离散格式,并验证静水平衡条件,确保控制方程各项离散结果的正确性;
步骤八、对经过上述步骤得到的控制方程各项进行一阶欧拉格式时间推进,采用完全隐式方案离散摩擦源项,通过分析得出隐式格式的解以显式方式更新流量通量,并通过Courant–Friedrichs–Levy条件限制时间步长;
步骤九、通过对时间步的不断推进,完成城市地表水流的数值模拟,并对模拟结果进行提取、处理和分析,得到城市地表水演进过程中的关键信息,包括但不限于水深、流速、最大水深、淹没历时等数据。
进一步地,步骤一中,获取的研究区域数字高程信息、土地利用类型遥感影像信息以及不透水建筑物的几何轮廓信息,需要限制为同一坐标系。其中地理数字高程信息用于构建研究区域数字高程模型;土地利用类型遥感影像信息用于获得研究区域内不同位置的糙率值,为数值模拟提供网格单元的初始曼宁系数值;不透水建筑物的几何轮廓信息用于对研究区域内不透水建筑进行概化,以描述建筑物的不透水特性。
进一步地,步骤二中,在用三角形非结构网格对研究区域进行空间离散时,将不透水建筑物的轮廓线作为网格剖分时的边界线,不透水建筑物所在区域不进行网格剖分,即不透水建筑物位置无网格并且不参与之后的数值计算;同时将得到的研究区域的数字高程信息及土地利用类型遥感影像信息在网格上插值得到每个网格处的高程值及糙率值。
进一步地,步骤三中,针对城市水流特点提出忽略动量项的简化型浅水方程组,并作为数值模拟的控制方程,分析原方程组各项,包括当地加速度项、压力项、对流项、底坡源项、摩擦源项代表的物理意义及作用,略去对流项,以简化数值模拟中大量的通量计算,得到简化型的二维浅水方程组:
式中,h为水深,u和v分别为x和y方向的流速,B为底坡高程,τ为摩擦项,g为重力加速度,τbx和τby分别为x和y方向的底坡摩擦力。
同时,将简化型浅水方程组的写成守恒格式如下:
式中U代表守恒变量,F和G分别为x和y方向的通量,Sb和Sf分别为底坡源项和摩擦源项,由曼宁公式得到各项具体形式如下:
其中n为曼宁系数。
进一步地,步骤四中,确定各物理变量,包括水深、流速、动量和质量在网格中的位置,采用单元中心式的非结构网格,所述简化型浅水方程组的离散后的积分形式如下:
其中,t为时间,下标i表示控制体的索引号,Ai为控制体Ωi的面积,n为单位外法向量,E=(F,G)T。
进一步,步骤五中,对求解简化型浅水方程组流量通量时遇到的黎曼问题,采用Roe格式黎曼求解器,构造相应的近似雅克比矩阵以及由其特征向量组成的相关系数矩阵,包括如下矩阵:
同时获得单元界面处各物理变量的黎曼状态,求解得出最终的通量计算公式如下:
进一步地,步骤六中,通过一种新的水面重构方法,修正各物理变量黎曼状态的表达式;水面重构方法的实现方式为,首先在单元界面处重建界面左右两侧的水面高程ηL和ηR:
ηL=ηP+max[0,min(bQ-bP,ηQ-ηP)];
ηR=ηQ+max[0,min(bP-bQ,ηP-ηQ)];
其中,η代表水面高程,b代表床底高程,bP为P单元的床底高程,bQ为Q单元的床底高程;ηP为P单元的水面高程,ηQ为Q单元的水面高程;
其次根据得到的水面高程重建单元界面左右两侧的床底高程bL和bR,并取左右侧的床底高程的最大值作为单元界面处床底高程bf:
bL=ηL-hP;
bR=ηR-hQ;
bf=max(bL,bR);
其中,hP为P单元处的水深,hQ为Q单元处的水深;因此单元界面处的黎曼状态量表示为:
ηL=max(0,ηL-bf);
hR=max(0,ηR-bf);
[hu]L=hLuP,[hv]L=hLvP;
[hu]R=hRuQ,[hv]R=hRvQ;
其中,hL为单元界面处左侧水深、hR为单元界面处右侧水深,[hu]L、[hv]L分别为单元界面处左侧的动量在x和y方向的分量,[hu]R、[hv]R分别为单元界面处右侧的动量在x和y方向的分量,代表P和Q单元中心x方向上的流速;代表P和Q单元中心y方向上的流速;同时对底坡源项在非结构网格上离散后,将得到的各黎曼状态值应用其中,因此底坡源项的迎风离散格式如下:
进一步地,步骤七中,通过修正单元界面处的高程值,处理城市地表水流数值模拟过程中的干湿水深动边问题,对bf进行如下修正并应用于所有计算单元:
其中:
式中εh为判定干单元的最小水深值,得到修正后的底坡源项的迎风离散格式如下:
进一步地,步骤八中,采用完全隐式方案离散摩擦源项,得到用于更新守恒物理量的显式计算方式:
进一步地,步骤八中,时间步长Δt严格受Courant–Friedrichs–Levy(CFL)条件限制:
其中,di是单元中心距离单元边缘的最小距离,CFL为0~1之间的非0值。ui为本单元中心位置的流速。
进一步地,步骤九中,完成城市地表水流的数值模拟后,得到数值模拟结果,并可根据相关需求(如预报需求、避险需求等),对结果进行进一步地分析和处理。
本发明的优点和积极效果在于:
本发明充分考虑城区洪水运动的特点后,简化现有数值模型的控制方程,大大降低了单元间数值通量部分的计算量,提高城市地表水流数值模拟的计算效率;在处理底坡源项时通过水面重构方法,解决了干湿水深变化产生的动边界问题,简化了处理干湿水深动边界问题的计算流程;使用隐式格式离散源项,并显式更新各守恒物理量,时间步长仅仅由CFL值决定,无需进一步约束。本发明提出的方法应用于城市地表水流数值模拟中,有效弥补现有数值模型模拟效率的不足,为城市地表水流数值模拟提供了一种新的途径。
附图说明
图1为本发明一种基于简化型浅水方程组的城市地表水流数值模拟方法的流程图;
图2为本发明实施例中一种基于简化型浅水方程组的城市地表水流数值模拟方法的在非结构网格上数值离散的示意图;
图3为本发明实施例中一种基于简化型浅水方程组的城市地表水流数值模拟方法的干湿水深动边界处理的示意图。
具体实施方式
以下结合附图及实施例对本发明的具体实施作进一步说明。下面参考的附图仅是示例性质的,本发明的实施方式不限于此。
如图1所示,一种基于简化型浅水方程组的城市地表水流数值模拟方法,包括以下步骤:
步骤一、获取研究区域的地理数字高程信息、土地利用类型遥感影像信息以及不透水建筑物的几何轮廓信息。这三者信息需要限制为同一坐标系,以保证三者在空间上联系的准确性,其中地理数字高程信息用于构建研究区域数字高程模型;土地利用类型遥感影像信息用于获得研究区域内不同位置的糙率值,为数值模拟提供网格单元的初始曼宁系数值;不透水建筑物的几何轮廓信息用于对研究区域内不透水建筑进行概化,以描述建筑物的不透水特性。
步骤二、采用三角形非结构网格对研究区域进行空间离散,将不透水建筑物的轮廓线作为网格剖分时的边界线,不透水建筑物所在区域不进行网格剖分,即不透水建筑物位置处无网格并且不参与之后的数值计算,使建筑物的不透水性在计算中得到体现;同时,将得到的研究区域的数字高程信息在网格上插值得到研究区域数字高程模型,每个网格的高程为数值模拟提供初始底坡高程值;将获取的土地利用类型遥感影像在网格上插值得到研究区域内不同位置的糙率值,为数值模拟提供每个网格单元的初始曼宁系数值。
步骤三、对完整的二维浅水方程组进行简化:针对城市水流特点提出忽略动量项的简化型浅水方程组,并作为数值模拟的控制方程,分析原方程组各项,包括当地加速度项、压力项、对流项、底坡源项、摩擦源项代表的物理意义及作用,略去对流项,以简化数值模拟中大量的通量计算,得到简化型的二维浅水方程组:
首先完整的浅水方程组如下:
其中第一个方程为质量守恒方程,第二个和第三个方程为动量守恒方程。式中,h为水深,u和v分别为x和y方向的流速,B为底坡高程,τ为摩擦项,g为重力加速度,τbx和τby分别为x和y方向的底坡摩擦力。方程组第二、三方程中第一项为局部加速度项,表示动量在任何固定位置的时间变化率,象征流体的不稳定性;第二、三项为对流加速度项,表示动量的空间梯度效应;第四、五项分别为底坡源项,表示流体运动过程中的重力作用;最后一项为摩擦源项。本发明提出对此方程组进行简化,忽略其中的对流加速度项,得到简化型浅水方程组如下:
接着,推导简化型浅水方程组的守恒格式及非线性矩阵形式,得到用于构建数值模型的控制方程。将简化型浅水方程写成矢量形式如下:
式中,t为时间,U代表守恒变量,F和G分别为x和y方向的通量,Sb和Sf分别为底坡源项和摩擦源项,由曼宁公式得到各项具体形式如下:
其中n为曼宁系数。
步骤四、为了适用城市地表特征,采用有限体积法对守恒格式简化型浅水方程组的控制方程在非结构网格上进行离散,最终得到简化型浅水方程组的离散后的积分格式。本实施例中,以控制体Ωi为例,首先得到控制方程的积分形式:
其中E=(F,G)T,应用高斯定理,上式变成:
式中,n为控制体的单位外法向量。本实施例中,积分的空间域被一组三角形单元所覆盖,但不一定与坐标方向对齐。在控制体Ωi中,应用上述公式进行离散逼近,其中,体积积分代表整个单元区域的积分,而面积分代表通过每个单元边界的总通量。用Ui表示当前时间步控制体Ωi上的守恒变量的平均值,对于每个单元得到:
式中,下标i表示控制体的索引号,Ai为控制体Ωi的面积。
同时,采用单元中心式的非结构网格,同时确定各物理变量,包括水深、流速、动量和质量在网格中的位置:本实施例中,当用有限体积法模拟水流流动时,首先要对计算区域进行网格划分,然后在选定的控制体内进行空间离散。如图2所示,P、Q为两个相邻单元的质心,阴影部分为单元P的控制体Ωi,w为两相邻单元的共有边。对于单元中心式网格,每个网格单元即为单一的控制体,网格单元互不重叠。每个控制体积上,针对控制方程确定了以单元为中心的有限体积方法,控制方程中的物理变量分布在网格单元内质心上,并且表示为分段常数。
步骤五、针对简化型浅水方程模型在数值通量计算中存在的黎曼问题,采用Roe格式黎曼求解器,通过Roe平均雅可比矩阵,解决单元界面处左右物理量不一致的问题。在步骤四中的简化型浅水方程组的离散后的积分形式中,第二项代表通过单元界面处的数值通量,对于每个网格计算单元边界的数值通量,通过边界中点法则逼近闭合曲线积分,可得:
其中k代表每个网格的边索引值,NE是每个单元的总边数,三角形网格单元中,NE=3,nk为第k条边的单位外法向量,lk为第k条边的长度。
为了评估单元P与相邻单元Q界面处相应的通量,在相邻网格单元共有边w左右两侧假设存在一维黎曼问题,对于单元P,假设在边界w左侧和右侧的数值通量分别为UL和UR。本实施例中,基于Roe格式黎曼求解器,首先构造恒定的雅克比矩阵,使得在单元边界处守恒物理量的差ΔU=UR-UL和通量差ΔE=ER-EL之间拥有一个满足守恒条件的线性关系:
其中第一条性质针对的是满足Rankine-Hugoniot条件的间断解;第二条性质是一种光滑性要求;第三条性质是为了保证问题的双曲性不变。由此,可以得到:
令P、Q两单元界面处的单位流量通量为(Ek·nk)=φPQ,可以得到单元边界处流量通量的计算方式为:
为了方便计算,流量通量可以写成如下形式:
步骤六、通过一种新的水面重构方法,重建黎曼状态量的表达式,并对简化型浅水方程组的底坡源项在网格上进行空间离散,得到底坡源项的迎风离散格式。同时由于Roe格式黎曼求解器不会自动满足静水平衡条件,因此通过对底坡源项的处理,可以使控制方程满足静水平衡条件,这一过程将在下一步实施步骤中进行验证。本实施例中,采用一种新的水面重构方法,对单元界面处两侧的物理变量进行重新定义,获得新的黎曼状态,克服了Roe格式黎曼求解器的局限性,满足静水平衡条件。水面重构方法的具体实施方式为:对于两个相邻单元P和Q,首先在单元界面处重建界面左右两侧的水面高程:
ηL=ηP+max[0,min(bQ-bP,ηQ-ηP)];
ηR=ηQ+max[0,min(bP-bQ,ηP-ηQ)];
式中,η代表水面高程,b代表床底高程。通过这一步的水面高程重建,我们可以得到界面处左右两侧的床底高程分别为:
bL=ηL-hP;
bR=ηR-hQ;
因此,单元界面处床底高程的黎曼状态量可以表示为:
bf=max(bL,bR);
同时,单元界面处左右两侧水深的黎曼状态量可以表示为:
hL=max(0,ηL-bf);
hR=max(0,ηR-bf);
此时,水深具有严格的非负性。随后流量通量计算公式中各黎曼状态量都得到了重建:
[hu]L=hLuP,[hv]L=hLvP;
[hu]R=hRuQ,[hv]R=hRvQ;
经过上述的水面重构方法,本实施例中考虑到了在表面流模拟过程中所有可能遇到的水位和床底高程的情况,并且在水面重构过程中,有效保留了单元内真正的床底变化情况,因此底坡源项可以使用下式简单的离散化:
步骤七、通过修正单元界面处高程值,处理城市地表水流数值模拟中的干湿水深动边界问题,得到修正后的底坡源项的迎风离散格式,并验证静水平衡条件。考虑到任意一个单元P,在模拟过程中,当其某个界面处的最终重建水面高程低于与该界面相邻单元界面处的床面高程时,可能会存在当水流到达一个干网格附近并无法继续演进时的情况。为了在干湿界面处保持静水平衡条件,以及床底高度突然变化引起的数值不稳定,对bf进行如下修正并应用于所有计算单元:
其中:
如图3所示,代表干湿边界交界处的情况,单元P为湿单元,Q为干单元,同时P单元的水位值低于Q单元的床底高程,在此情况下Δb=bf-ηP,经过水面重构后即会降低界面处右侧的床底高程值至左侧水位值,这样可以有效地保持干湿界面处的静水平衡。下面证明本发明采用的水面重构法可以满足静水平衡条件,对于一个静止水面可得:
其中前两项经过代数运算后结果为0,表达式变为:
这与静水条件下的数值通量计算结果相一致,则净水平衡条件得到满足,同时确保了控制方程各项离散结果的正确性。
步骤八、对经过上述步骤得到的控制方程各项进行一阶欧拉格式时间推进,采用完全隐式方案离散摩擦源项,通过分析得出隐式格式的解以显式方式更新流量通量,并通过Courant–Friedrichs–Levy条件限制时间步长,完成对城市地表水流数值的模拟。本实施例中,控制方程的一阶欧拉时间推进公式如下:
其中上标j表示当前时间步,Δt为时间步长,Ai为网格单元面积,i表示控制体的索引号。数值方案中,流量通量和底坡源项基于当前时间步j离散,而摩擦源项基于下一时间步隐式离散,但是实际计算时,通过推导有效的显示公式来实现。为了推导出显示的表达式,对时间推进公式改写成如下形式:
将两式相除,得到:
本实施例中,以x方向上分量为例,带入x方向上时间推进公式,得到:
同时为了避免hj为0引起的非法计算,修正后得到:
经过上述推导,摩擦源项的隐式离散格式,可以被显式地计算并应用于简化型浅水方程组。本发明提出的方法模型稳定模拟的时间步长仅受CFL条件限制:
其中,di是单元中心距离单元边缘的最小距离,CFL为0~1之间的非0值。至此,简化型浅水方程组中的各项在空间上都已经完成离散。并且在时间上也完成了向下一时间步的推进,可直接进行数值模拟计算。
步骤九、通过对时间步的不断推进,完成城市地表水流的数值模拟,并对模拟结果进行提取、处理和分析,可以得到城市地表水演进过程中的关键信息,包括但不限于水深、流速、最大水深、淹没历时等数据。
在本说明书的叙述中提到的相同或类似的符号和标注,代表相同或近似的物理意义或具有相同或近似的功能,并且本说明书中使用的图例,仅仅是为了更好的解释本发明,本发明的适用性并不限制于此。
本发明一种基于简化型浅水方程组的城市地表水流数值模拟方法,根据说明书提供的示例及相关步骤,可以在城市地区实现复杂水流运动的数值模拟,简化传统城市水流模型的相关模拟流程,提高城市地表水流数值模拟的效率,为构建城市洪涝预警预报系统的数值模型部分提供一种新的途径,实现对城市洪水险情的有效预警以及规避损失。
Claims (3)
1.一种基于简化型浅水方程组的城市地表水流数值模拟方法,其特征在于,包括以下步骤:
步骤一、获取研究区域的地理数字高程信息、土地利用类型遥感影像信息以及不透水建筑物的几何轮廓信息;
步骤二、采用三角形非结构网格对研究区域进行空间离散,其中不透水建筑物部分不进行网格剖分,将步骤一获得的研究区域的数字高程信息及土地利用类型遥感影像信息在网格上插值得到每个网格的高程值及糙率值;
步骤三、对完整的二维浅水方程组进行简化,针对城市水流特点推导出一种忽略动量方程中对流加速度项的简化型浅水方程组的守恒格式及非线性矩阵形式,将简化型浅水方程组作为用于构建数值模型的控制方程,并基于所述控制方程对城市地表水流进行数值模拟;
通过分析完整的二维浅水方程组各项,包括当地加速度项、压力项、对流项、底坡源项、摩擦源项代表的物理意义及作用,略去对流项,以简化数值模拟中大量的通量计算,得到的忽略动量方程中对流加速度项的简化型浅水方程组如下:
式中,h为水深,u和v分别为x和y方向的流速,t为时间,B为底坡高程,τ为摩擦项,g为重力加速度,τbx和τby分别为x和y方向的底坡摩擦力;
简化型浅水方程组的守恒格式如下:
式中U代表守恒变量,F和G分别为x和y方向的通量,Sb和Sf分别为底坡源项和摩擦源项,由曼宁公式得到各项具体形式如下:
其中n为曼宁系数;
步骤四、为了适用城市地表特征,采用有限体积法对守恒格式简化型浅水方程组的控制方程组在三角形非结构网格上进行离散,同时确定各物理变量,包括水深、流速、动量和质量在网格中的位置,最终得到简化型浅水方程组的离散后的积分格式,具体如下:
步骤五、针对简化型浅水方程数值模型单元间的流量通量计算中存在的黎曼问题,采用Roe格式黎曼求解器,构造相应的近似雅克比矩阵以及由其特征值和特征向量组成的相关系数矩阵,包括如下矩阵:
获得单元界面处各物理变量的黎曼状态,以解决单元界面处左右物理量不一致的问题,求解得出最终的通量计算公式如下:
式中,φPQ为单元间的流量通量,hL和hR分别表示单元界面左右两侧的水深,Δh=hR-hL,uL和uR分别表示单元界面处左右两侧在x方向上的流速,vL和vR分别表示单元界面左右两侧在y方向上的流速;
步骤六、通过一种新的水面重构方法,修正黎曼状态量的表达式,面重构方法的实现方式为,首先在单元界面处重建界面左右两侧的水面高程ηL和ηR:
ηL=ηP+max[0,min(bQ-bP,ηQ-ηP)];
ηR=ηQ+max[0,min(bP-bQ,ηP-ηQ)];
其中,η代表水面高程,b代表床底高程,bP为P单元的床底高程,bQ为Q单元的床底高程;ηP为P单元的水面高程,ηQ为Q单元的水面高程;
其次根据得到的水面高程重建单元界面左右两侧的床底高程bL和bR,并取左右侧的床底高程的最大值作为单元界面处床底高程bf:
bL=ηL-hP;
bR=ηR-hQ;
bf=max(bL,bR);
其中,hP为P单元处的水深,hQ为Q单元处的水深,bf为单元界面处的床底高程值;因此单元界面处的黎曼状态量表示为:
hL=max(0,ηL-bf);
hR=max(0,ηR-bf);
[hu]L=hLuP,[hv]L=hLvP;
[hu]R=hRuQ,[hv]R=hRvQ;
其中,[hu]L、[hv]L分别为单元界面处左侧的动量在x和y方向的分量,[hu]R、[hv]R分别为单元界面处右侧的动量在x和y方向的分量,代表P和Q单元中心x方向上的流速,[hu]P和[hu]Q分别代表P和Q单元中心x方向上的动量;代表P和Q单元中心y方向上的流速,[hv]P和[hv]Q分别代表P和Q单元中心y方向上的动量;
并对简化型浅水方程组的底坡源项在非结构网格上离散后,将得到的各黎曼状态值应用其中,因此底坡源项的迎风离散格式如下:
步骤七、通过修正单元界面处的高程值,处理城市地表水流数值模拟过程中的干湿水深动边问题,对bf进行如下修正并应用于所有计算单元:
其中:
步骤八、对经过上述步骤得到的控制方程各项进行一阶欧拉格式时间推进,采用完全隐式方案离散摩擦源项,得到用于更新守恒物理量的显式计算方式:
式中,x和y表示物理量在x和y方向上的分量,有着上标j的物理量代表物理量在当前时间步的值,和分别表示下一时间步内单元物理变量在x和y方向上的分量,hj为当前时间步内单元水深,Δt为时间步长,其中,Ax和Ay分别表示动量在x和y方向上的分量,表示当前时间步内单元第k条边界面处的通量,和分别表示当前时间步内单元物理变量在x和y方向上的分量,Sj bi表示当前时间步内单元i的底坡源项;
时间步长Δt严格受Courant-Friedrichs-Levy(CFL)条件限制:
其中,di是单元中心距离单元边缘的最小距离,CFL为0~1之间的非0值,ui为本单元中心位置的流速;
步骤九、通过对时间步的不断推进,完成城市地表水流的数值模拟。
2.根据权利要求1所述的一种基于简化型浅水方程组的城市地表水流数值模拟方法,其特征在于,步骤九还包括对模拟结果进行提取、处理和分析,得到城市地表水演进过程中关键信息,关键信息包括水深、流速、最大水深、淹没历时数据。
3.根据权利要求1所述的一种基于简化型浅水方程组的城市地表水流数值模拟方法,其特征在于,步骤一中,获取的研究区域数字高程信息、土地利用类型遥感影像信息以及不透水建筑物的几何轮廓信息,需要限制为同一坐标系;
步骤二中,在用三角形非结构网格对研究区域进行空间离散时,将不透水建筑物的轮廓线作为网格剖分时的边界线,不透水建筑物所在区域不进行网格剖分,即不透水建筑物位置无网格并且不参与之后的数值计算;同时将得到的研究区域的数字高程信息及土地利用类型遥感影像信息在网格上插值得到每个网格处的高程值及糙率值。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201911147817.0A CN110765694B (zh) | 2019-11-21 | 2019-11-21 | 一种基于简化型浅水方程组的城市地表水流数值模拟方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201911147817.0A CN110765694B (zh) | 2019-11-21 | 2019-11-21 | 一种基于简化型浅水方程组的城市地表水流数值模拟方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110765694A CN110765694A (zh) | 2020-02-07 |
CN110765694B true CN110765694B (zh) | 2021-11-05 |
Family
ID=69338854
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201911147817.0A Active CN110765694B (zh) | 2019-11-21 | 2019-11-21 | 一种基于简化型浅水方程组的城市地表水流数值模拟方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110765694B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111931263B (zh) * | 2020-05-19 | 2021-09-14 | 长江水利委员会水文局 | 一种基于优化求解平面二维浅水方程的恒定流模拟方法 |
CN112257313B (zh) * | 2020-10-21 | 2024-05-14 | 西安理工大学 | 一种基于gpu加速的污染物输移高分辨率数值模拟方法 |
CN115130264B (zh) * | 2022-09-01 | 2022-11-25 | 浙江远算科技有限公司 | 一种基于径流耦合仿真的城市内涝预测方法及系统 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107784141A (zh) * | 2016-08-31 | 2018-03-09 | 施勇 | 一种二维有限控制体积计算的加速方法 |
CN108446502A (zh) * | 2018-03-22 | 2018-08-24 | 中国水利水电科学研究院 | 一种利用完整二维浅水方程组获得流域单位线的方法 |
CN109492259A (zh) * | 2018-10-15 | 2019-03-19 | 华北水利水电大学 | 一种城市水文模拟系统 |
CN109543275A (zh) * | 2018-11-15 | 2019-03-29 | 中国水利水电科学研究院 | 一种城区地表径流二维数值模拟方法 |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7239990B2 (en) * | 2003-02-20 | 2007-07-03 | Robert Struijs | Method for the numerical simulation of a physical phenomenon with a preferential direction |
CN106940840A (zh) * | 2017-03-14 | 2017-07-11 | 东南大学 | 一种城市内涝灾害风险评估方法 |
CN108920799B (zh) * | 2018-06-22 | 2022-04-26 | 中国科学院地理科学与资源研究所 | 基于正方形自适应网格的二维水文-水动力单向耦合方法 |
-
2019
- 2019-11-21 CN CN201911147817.0A patent/CN110765694B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107784141A (zh) * | 2016-08-31 | 2018-03-09 | 施勇 | 一种二维有限控制体积计算的加速方法 |
CN108446502A (zh) * | 2018-03-22 | 2018-08-24 | 中国水利水电科学研究院 | 一种利用完整二维浅水方程组获得流域单位线的方法 |
CN109492259A (zh) * | 2018-10-15 | 2019-03-19 | 华北水利水电大学 | 一种城市水文模拟系统 |
CN109543275A (zh) * | 2018-11-15 | 2019-03-29 | 中国水利水电科学研究院 | 一种城区地表径流二维数值模拟方法 |
Non-Patent Citations (1)
Title |
---|
Haijun Yu, et al..Efficient Finite-Volume Model for Shallow-Water Flows Using an Implicit Dual Time-Stepping Method.《Journal of Hydraulic Engineering》.2015, * |
Also Published As
Publication number | Publication date |
---|---|
CN110765694A (zh) | 2020-02-07 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110765694B (zh) | 一种基于简化型浅水方程组的城市地表水流数值模拟方法 | |
JP5209298B2 (ja) | 流れのシミュレーション計算方法およびシステム | |
US10439594B2 (en) | Actually-measured marine environment data assimilation method based on sequence recursive filtering three-dimensional variation | |
Mandli et al. | Adaptive mesh refinement for storm surge | |
CN110232471B (zh) | 一种降水传感网节点布局优化方法及装置 | |
CN111767597B (zh) | 一种城市模型验证方法、装置、设备及存储介质 | |
Yue et al. | An adaptive method of high accuracy surface modeling and its application to simulating elevation surfaces | |
CN109726465B (zh) | 基于非结构曲边网格的三维无粘低速绕流的数值模拟方法 | |
CN107133686A (zh) | 基于时空数据模型的城市级pm2.5浓度预测方法 | |
Ferragut et al. | Comparison between 2.5-D and 3-D realistic models for wind field adjustment | |
CN109726433A (zh) | 基于曲面边界条件的三维无粘低速绕流的数值模拟方法 | |
CN104749625B (zh) | 一种基于正则化技术的地震数据倾角估计方法及装置 | |
CN110990926B (zh) | 一种基于面积修正率的城市地表建筑水动力学仿真方法 | |
Han et al. | CoolVox: Advanced 3D convolutional neural network models for predicting solar radiation on building facades | |
CN109033181B (zh) | 一种复杂地形地区风场地理数值模拟方法 | |
CN117951989A (zh) | 一种降雨边坡临界滑面搜索方法、装置及可读存储介质 | |
Tymków et al. | 3D GIS for flood modelling in river valleys | |
Li et al. | An orthogonal terrain-following coordinate and its preliminary tests using 2-D idealized advection experiments | |
WO2019128018A1 (zh) | 确定风资源的方法和装置 | |
US9633165B1 (en) | Composite design direction | |
CN111898819B (zh) | 空间网格划分方法及装置 | |
KR102026153B1 (ko) | 행위자 기반 모형을 포함하는 다중 규모 모형을 이용한 남조류 이송·확산 거동의 수치모의 방법 | |
CN113033035A (zh) | 一种污染物扩散区域动态模拟方法、系统及装置 | |
CN114004175B (zh) | 一种快速查找全域壁面距离和无量纲壁面距离的方法 | |
Liu et al. | 3D seabed terrain establishment based on moving fractal interpolation |
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 |