CN110765694B - 一种基于简化型浅水方程组的城市地表水流数值模拟方法 - Google Patents

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方向的通量，S_b和S_f分别为底坡源项和摩擦源项，由曼宁公式得到各项具体形式如下：

其中n为曼宁系数。

进一步地，步骤四中，确定各物理变量，包括水深、流速、动量和质量在网格中的位置，采用单元中心式的非结构网格，所述简化型浅水方程组的离散后的积分形式如下：

其中，t为时间，下标i表示控制体的索引号，A_i为控制体Ω_i的面积，n为单位外法向量，E＝(F,G)^T。

进一步，步骤五中，对求解简化型浅水方程组流量通量时遇到的黎曼问题，采用Roe格式黎曼求解器，构造相应的近似雅克比矩阵以及由其特征向量组成的相关系数矩阵，包括如下矩阵：

雅克比矩阵的近似Roe平均矩阵

矩阵

的特征向量矩阵

及其逆矩阵：

的改进形式

式中，下标x和y分别表示物理量在x和y方向的分量，

为Roe平均水深；

同时获得单元界面处各物理变量的黎曼状态，求解得出最终的通量计算公式如下：

式中，φ_PQ为单元间的流量通量，

Δh＝h_R-h_L，下标L和R分别表示单元界面左右侧。

进一步地，步骤六中，通过一种新的水面重构方法，修正各物理变量黎曼状态的表达式；水面重构方法的实现方式为，首先在单元界面处重建界面左右两侧的水面高程η_L和η_R：

η_L＝η_P+max[0,min(b_Q-b_P,η_Q-η_P)]；

η_R＝η_Q+max[0,min(b_P-b_Q,η_P-η_Q)]；

其中，η代表水面高程，b代表床底高程，b_P为P单元的床底高程，b_Q为Q单元的床底高程；η_P为P单元的水面高程，η_Q为Q单元的水面高程；

其次根据得到的水面高程重建单元界面左右两侧的床底高程b_L和b_R，并取左右侧的床底高程的最大值作为单元界面处床底高程b_f：

b_L＝η_L-h_P；

b_R＝η_R-h_Q；

b_f＝max(b_L,b_R)；

其中，h_P为P单元处的水深，h_Q为Q单元处的水深；因此单元界面处的黎曼状态量表示为：

η_L＝max(0,η_L-b_f)；

h_R＝max(0,η_R-b_f)；

[hu]_L＝h_Lu_P,[hv]_L＝h_Lv_P；

[hu]_R＝h_Ru_Q,[hv]_R＝h_Rv_Q；

其中，h_L为单元界面处左侧水深、h_R为单元界面处右侧水深，[hu]_L、[hv]_L分别为单元界面处左侧的动量在x和y方向的分量，[hu]_R、[hv]_R分别为单元界面处右侧的动量在x和y方向的分量，

代表P和Q单元中心x方向上的流速；

代表P和Q单元中心y方向上的流速；同时对底坡源项在非结构网格上离散后，将得到的各黎曼状态值应用其中，因此底坡源项的迎风离散格式如下：

式中，

为底坡源项，b_P为P单元的床底高程，k为单元边的索引号，h_L,k即为P单元在第k条边界面处左侧的水深，n_k和l_k分别是第k条边的单位外法向量及长度，b_f为单元界面处的床底高程值。

进一步地，步骤七中，通过修正单元界面处的高程值，处理城市地表水流数值模拟过程中的干湿水深动边问题，对b_f进行如下修正并应用于所有计算单元：

其中：

式中ε_h为判定干单元的最小水深值，得到修正后的底坡源项的迎风离散格式如下：

式中，

为修正后界面处的高程值。

进一步地，步骤八中，采用完全隐式方案离散摩擦源项，得到用于更新守恒物理量的显式计算方式：

式中，x和y表示物理量在x和y方向上的分量，有着上标j的物理量代表物理量在当前时间步的值，Δt为时间步长，其中m＝U^j+ΔtA，

进一步地，步骤八中，时间步长Δt严格受Courant–Friedrichs–Levy(CFL)条件限制：

其中，d_i是单元中心距离单元边缘的最小距离，CFL为0～1之间的非0值。u_i为本单元中心位置的流速。

进一步地，步骤九中，完成城市地表水流的数值模拟后，得到数值模拟结果，并可根据相关需求(如预报需求、避险需求等)，对结果进行进一步地分析和处理。

本发明的优点和积极效果在于：

本发明充分考虑城区洪水运动的特点后，简化现有数值模型的控制方程，大大降低了单元间数值通量部分的计算量，提高城市地表水流数值模拟的计算效率；在处理底坡源项时通过水面重构方法，解决了干湿水深变化产生的动边界问题，简化了处理干湿水深动边界问题的计算流程；使用隐式格式离散源项，并显式更新各守恒物理量，时间步长仅仅由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中，应用上述公式进行离散逼近，其中，体积积分代表整个单元区域的积分，而面积分代表通过每个单元边界的总通量。用U_i表示当前时间步控制体Ω_i上的守恒变量的平均值，对于每个单元得到：

式中，下标i表示控制体的索引号，A_i为控制体Ω_i的面积。

同时，采用单元中心式的非结构网格，同时确定各物理变量，包括水深、流速、动量和质量在网格中的位置：本实施例中，当用有限体积法模拟水流流动时，首先要对计算区域进行网格划分，然后在选定的控制体内进行空间离散。如图2所示，P、Q为两个相邻单元的质心，阴影部分为单元P的控制体Ω_i，w为两相邻单元的共有边。对于单元中心式网格，每个网格单元即为单一的控制体，网格单元互不重叠。每个控制体积上，针对控制方程确定了以单元为中心的有限体积方法，控制方程中的物理变量分布在网格单元内质心上，并且表示为分段常数。

步骤五、针对简化型浅水方程模型在数值通量计算中存在的黎曼问题，采用Roe格式黎曼求解器，通过Roe平均雅可比矩阵，解决单元界面处左右物理量不一致的问题。在步骤四中的简化型浅水方程组的离散后的积分形式中，第二项代表通过单元界面处的数值通量，对于每个网格计算单元边界的数值通量，通过边界中点法则逼近闭合曲线积分，可得：

其中k代表每个网格的边索引值，NE是每个单元的总边数，三角形网格单元中，NE＝3，n_k为第k条边的单位外法向量，l_k为第k条边的长度。

为了评估单元P与相邻单元Q界面处相应的通量，在相邻网格单元共有边w左右两侧假设存在一维黎曼问题，对于单元P，假设在边界w左侧和右侧的数值通量分别为U_L和U_R。本实施例中，基于Roe格式黎曼求解器，首先构造恒定的雅克比矩阵，使得在单元边界处守恒物理量的差ΔU＝U_R-U_L和通量差ΔE＝E_R-E_L之间拥有一个满足守恒条件的线性关系：

其中

为常数矩阵，称为雅可比矩阵J的Roe平均矩阵，它具有U性质：

(1)

(2)

(3)

有实数特征根以及完备的线性无关的特征向量组。

其中第一条性质针对的是满足Rankine-Hugoniot条件的间断解；第二条性质是一种光滑性要求；第三条性质是为了保证问题的双曲性不变。由此，可以得到：

根据Roe格式黎曼求解器的要求，矩阵

和J有相同的形状，矩阵中的值由

及

确定，于是得到：

其中上标～表示物理量的Roe平均值，

同时可以得到

矩阵的特征向量和特征值：

令P、Q两单元界面处的单位流量通量为(E_k·n_k)＝φ_PQ，可以得到单元边界处流量通量的计算方式为：

其中

是由特征值组成的绝对值对角矩阵，并且

为了方便计算，流量通量可以写成如下形式：

其中

得到:

将

代入通量计算公式，最终得到单元P中与单元Q相邻界面处的流量通量为：

步骤六、通过一种新的水面重构方法,重建黎曼状态量的表达式，并对简化型浅水方程组的底坡源项在网格上进行空间离散，得到底坡源项的迎风离散格式。同时由于Roe格式黎曼求解器不会自动满足静水平衡条件，因此通过对底坡源项的处理，可以使控制方程满足静水平衡条件，这一过程将在下一步实施步骤中进行验证。本实施例中，采用一种新的水面重构方法，对单元界面处两侧的物理变量进行重新定义，获得新的黎曼状态，克服了Roe格式黎曼求解器的局限性，满足静水平衡条件。水面重构方法的具体实施方式为：对于两个相邻单元P和Q，首先在单元界面处重建界面左右两侧的水面高程：

η_L＝η_P+max[0,min(b_Q-b_P,η_Q-η_P)]；

η_R＝η_Q+max[0,min(b_P-b_Q,η_P-η_Q)]；

式中，η代表水面高程，b代表床底高程。通过这一步的水面高程重建，我们可以得到界面处左右两侧的床底高程分别为：

b_L＝η_L-h_P；

b_R＝η_R-h_Q；

因此，单元界面处床底高程的黎曼状态量可以表示为：

b_f＝max(b_L,b_R)；

同时，单元界面处左右两侧水深的黎曼状态量可以表示为：

h_L＝max(0,η_L-b_f)；

h_R＝max(0,η_R-b_f)；

此时，水深具有严格的非负性。随后流量通量计算公式中各黎曼状态量都得到了重建：

[hu]_L＝h_Lu_P,[hv]_L＝h_Lv_P；

[hu]_R＝h_Ru_Q,[hv]_R＝h_Rv_Q；

其中，

代表P和Q网格单元中心x方向上的流速，y方向上的流速与其类似，本实施例中，重新定义了流量通量项中的各物理变量黎曼状态的表达式，并将其应用在计算流量通量项的离散格式中。

经过上述的水面重构方法，本实施例中考虑到了在表面流模拟过程中所有可能遇到的水位和床底高程的情况，并且在水面重构过程中，有效保留了单元内真正的床底变化情况，因此底坡源项可以使用下式简单的离散化：

式中，

为底坡源项，b_P为P单元的床底高程，A_i为网格单元面积，k为单元边的索引号，n_k和l_k分别是第k条边的单位外法向量及长度，b_f为单元界面处的床底高程值。

步骤七、通过修正单元界面处高程值，处理城市地表水流数值模拟中的干湿水深动边界问题，得到修正后的底坡源项的迎风离散格式，并验证静水平衡条件。考虑到任意一个单元P，在模拟过程中，当其某个界面处的最终重建水面高程低于与该界面相邻单元界面处的床面高程时，可能会存在当水流到达一个干网格附近并无法继续演进时的情况。为了在干湿界面处保持静水平衡条件，以及床底高度突然变化引起的数值不稳定，对b_f进行如下修正并应用于所有计算单元：

其中：

式中ε_h为判定干单元的最小水深值，

为修正后界面处的高程值。针对城市地表水流数值模拟中的干湿水深动边界问题进行修正后的底坡源项离散格式如下：

如图3所示，代表干湿边界交界处的情况，单元P为湿单元，Q为干单元，同时P单元的水位值低于Q单元的床底高程，在此情况下Δb＝b_f-η_P，经过水面重构后

即会降低界面处右侧的床底高程值至左侧水位值，这样可以有效地保持干湿界面处的静水平衡。下面证明本发明采用的水面重构法可以满足静水平衡条件，对于一个静止水面

可得：

的动量分量可以通过上述公式写成：

其中前两项经过代数运算后结果为0，表达式变为：

这与静水条件下的数值通量计算结果相一致，则净水平衡条件得到满足，同时确保了控制方程各项离散结果的正确性。

步骤八、对经过上述步骤得到的控制方程各项进行一阶欧拉格式时间推进，采用完全隐式方案离散摩擦源项，通过分析得出隐式格式的解以显式方式更新流量通量，并通过Courant–Friedrichs–Levy条件限制时间步长，完成对城市地表水流数值的模拟。本实施例中，控制方程的一阶欧拉时间推进公式如下：

其中上标j表示当前时间步，Δt为时间步长，A_i为网格单元面积，i表示控制体的索引号。数值方案中，流量通量和底坡源项基于当前时间步j离散，而摩擦源项基于下一时间步隐式离散，但是实际计算时，通过推导有效的显示公式来实现。为了推导出显示的表达式，对时间推进公式改写成如下形式：

式中下标x和y分别表示物理量在x和y方向上的分量，A_x和A_y分别表示动量

在x和y方向上的分量。采用下述方法得到上式中

和

的解，首先修改上述两式为：

将两式相除，得到：

本实施例中，以x方向上分量为例，带入x方向上时间推进公式，得到：

其中

值得注意的是根据

的正负值，上式中可以得到多达四个解，然而只有一个解是符合流动的物理特征。四个根中哪一个是符合物理规律的，取决于m_x和

根据对m_x和

的正负性相一致的特点，得到正确解的表达式为：

如果h^j的值很小，接近于干水深判定值，

的可能由于数值过大导致数值误差，为了避免这种情况，对其平方根处理，得到下式：

同时为了避免h^j为0引起的非法计算，修正后得到：

同理，对于

得到：

经过上述推导，摩擦源项的隐式离散格式，可以被显式地计算并应用于简化型浅水方程组。本发明提出的方法模型稳定模拟的时间步长仅受CFL条件限制：

其中，d_i是单元中心距离单元边缘的最小距离，CFL为0～1之间的非0值。至此，简化型浅水方程组中的各项在空间上都已经完成离散。并且在时间上也完成了向下一时间步的推进，可直接进行数值模拟计算。

步骤九、通过对时间步的不断推进，完成城市地表水流的数值模拟，并对模拟结果进行提取、处理和分析，可以得到城市地表水演进过程中的关键信息，包括但不限于水深、流速、最大水深、淹没历时等数据。

在本说明书的叙述中提到的相同或类似的符号和标注，代表相同或近似的物理意义或具有相同或近似的功能，并且本说明书中使用的图例，仅仅是为了更好的解释本发明，本发明的适用性并不限制于此。

本发明一种基于简化型浅水方程组的城市地表水流数值模拟方法，根据说明书提供的示例及相关步骤，可以在城市地区实现复杂水流运动的数值模拟，简化传统城市水流模型的相关模拟流程，提高城市地表水流数值模拟的效率，为构建城市洪涝预警预报系统的数值模型部分提供一种新的途径，实现对城市洪水险情的有效预警以及规避损失。

Claims (3)

1.一种基于简化型浅水方程组的城市地表水流数值模拟方法，其特征在于，包括以下步骤：

步骤一、获取研究区域的地理数字高程信息、土地利用类型遥感影像信息以及不透水建筑物的几何轮廓信息；

步骤二、采用三角形非结构网格对研究区域进行空间离散，其中不透水建筑物部分不进行网格剖分，将步骤一获得的研究区域的数字高程信息及土地利用类型遥感影像信息在网格上插值得到每个网格的高程值及糙率值；

步骤三、对完整的二维浅水方程组进行简化，针对城市水流特点推导出一种忽略动量方程中对流加速度项的简化型浅水方程组的守恒格式及非线性矩阵形式，将简化型浅水方程组作为用于构建数值模型的控制方程，并基于所述控制方程对城市地表水流进行数值模拟；

通过分析完整的二维浅水方程组各项，包括当地加速度项、压力项、对流项、底坡源项、摩擦源项代表的物理意义及作用，略去对流项，以简化数值模拟中大量的通量计算，得到的忽略动量方程中对流加速度项的简化型浅水方程组如下：

式中，h为水深，u和v分别为x和y方向的流速，t为时间，B为底坡高程，τ为摩擦项，g为重力加速度，τ_bx和τ_by分别为x和y方向的底坡摩擦力；

简化型浅水方程组的守恒格式如下：

式中U代表守恒变量，F和G分别为x和y方向的通量，S_b和S_f分别为底坡源项和摩擦源项，由曼宁公式得到各项具体形式如下：

其中n为曼宁系数；

步骤四、为了适用城市地表特征，采用有限体积法对守恒格式简化型浅水方程组的控制方程组在三角形非结构网格上进行离散，同时确定各物理变量，包括水深、流速、动量和质量在网格中的位置，最终得到简化型浅水方程组的离散后的积分格式，具体如下：

其中，下标i表示控制体的索引号，Ω代表控制体，A_i为控制体Ω_i的面积，U_i表示当前时间步控制体Ω_i上的守恒变量的平均值，

为单位外法向量，E＝(F，G)^T，S＝(S_b，S_f)^T；

步骤五、针对简化型浅水方程数值模型单元间的流量通量计算中存在的黎曼问题，采用Roe格式黎曼求解器，构造相应的近似雅克比矩阵以及由其特征值和特征向量组成的相关系数矩阵，包括如下矩阵：

雅克比矩阵的近似Roe平均矩阵

矩阵

的特征向量矩阵

及其逆矩阵：

的改进形式

式中，

为Roe平均水深，

h_L和h_R分别表示单元界面左右两侧的水深，

和

分别表示单位外法向量在x和y方向上的分量；

获得单元界面处各物理变量的黎曼状态，以解决单元界面处左右物理量不一致的问题，求解得出最终的通量计算公式如下：

式中，φ_PQ为单元间的流量通量，h_L和h_R分别表示单元界面左右两侧的水深，Δh＝h_R-h_L，u_L和u_R分别表示单元界面处左右两侧在x方向上的流速，v_L和v_R分别表示单元界面左右两侧在y方向上的流速；

步骤六、通过一种新的水面重构方法，修正黎曼状态量的表达式，面重构方法的实现方式为，首先在单元界面处重建界面左右两侧的水面高程η_L和η_R：

η_L＝η_P+max[0，min(b_Q-b_P，η_Q-η_P)]；

η_R＝η_Q+max[0，min(b_P-b_Q，η_P-η_Q)]；

其中，η代表水面高程，b代表床底高程，b_P为P单元的床底高程，b_Q为Q单元的床底高程；η_P为P单元的水面高程，η_Q为Q单元的水面高程；

其次根据得到的水面高程重建单元界面左右两侧的床底高程b_L和b_R，并取左右侧的床底高程的最大值作为单元界面处床底高程b_f：

b_L＝η_L-h_P；

b_R＝η_R-h_Q；

b_f＝max(b_L，b_R)；

其中，h_P为P单元处的水深，h_Q为Q单元处的水深，b_f为单元界面处的床底高程值；因此单元界面处的黎曼状态量表示为：

h_L＝max(0，η_L-b_f)；

h_R＝max(0，η_R-b_f)；

[hu]_L＝h_Lu_P，[hv]_L＝h_Lv_P；

[hu]_R＝h_Ru_Q，[hv]_R＝h_Rv_Q；

其中，[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方向上的动量；

并对简化型浅水方程组的底坡源项在非结构网格上离散后，将得到的各黎曼状态值应用其中，因此底坡源项的迎风离散格式如下：

式中，

为底坡源项，b_P为P单元的床底高程，k为单元边的索引号，h_L，k即为P单元在第k条边界面处左侧的水深，b_f，k即为P单元在第k条边界面处的床底高程，

和l_k分别是第k条边的单位外法向量及长度；

步骤七、通过修正单元界面处的高程值，处理城市地表水流数值模拟过程中的干湿水深动边问题，对b_f进行如下修正并应用于所有计算单元：

其中：

式中，

为修正后界面处的高程值，ε_h为判定干单元的最小水深值，得到修正后的底坡源项的迎风离散格式如下：

式中，

为修正后第k条边界面处的床底高程，并验证静水平衡条件，确保控制方程各项离散结果的正确性；

步骤八、对经过上述步骤得到的控制方程各项进行一阶欧拉格式时间推进，采用完全隐式方案离散摩擦源项，得到用于更新守恒物理量的显式计算方式：

式中，x和y表示物理量在x和y方向上的分量，有着上标j的物理量代表物理量在当前时间步的值，

和

分别表示下一时间步内单元物理变量在x和y方向上的分量，h^j为当前时间步内单元水深，Δt为时间步长，

其中，A_x和A_y分别表示动量

在x和y方向上的分量，

表示当前时间步内单元第k条边界面处的通量，

和

分别表示当前时间步内单元物理变量在x和y方向上的分量，S^j _bi表示当前时间步内单元i的底坡源项；

时间步长Δt严格受Courant-Friedrichs-Levy(CFL)条件限制：

其中，d_i是单元中心距离单元边缘的最小距离，CFL为0～1之间的非0值，u_i为本单元中心位置的流速；

步骤九、通过对时间步的不断推进，完成城市地表水流的数值模拟。

2.根据权利要求1所述的一种基于简化型浅水方程组的城市地表水流数值模拟方法，其特征在于，步骤九还包括对模拟结果进行提取、处理和分析，得到城市地表水演进过程中关键信息，关键信息包括水深、流速、最大水深、淹没历时数据。

3.根据权利要求1所述的一种基于简化型浅水方程组的城市地表水流数值模拟方法，其特征在于，步骤一中，获取的研究区域数字高程信息、土地利用类型遥感影像信息以及不透水建筑物的几何轮廓信息，需要限制为同一坐标系；

步骤二中，在用三角形非结构网格对研究区域进行空间离散时，将不透水建筑物的轮廓线作为网格剖分时的边界线，不透水建筑物所在区域不进行网格剖分，即不透水建筑物位置无网格并且不参与之后的数值计算；同时将得到的研究区域的数字高程信息及土地利用类型遥感影像信息在网格上插值得到每个网格处的高程值及糙率值。

Images