一种基于显卡加速耦合管网的高精度城市雨洪模拟方法

本发明属于数值模拟技术领域，涉及一种基于显卡加速耦合管网的高精度城市雨洪模拟方法。

背景技术：

洪涝灾害是在强降水、外江洪水、潮水顶托、蓄滞不足和排水不畅等因素作用下导致河湖水位上涨或雨水地表径流积聚而产生的地表积水现象。由于城市化的快速发展，可透水下垫面面积比例不断减小，近年来，全球气候变化、极端降雨等多因素影响下的暴雨洪涝灾害日益频发，已成为制约城市发展与运行管理的一个非常突出的问题。因此，研究城市雨洪的水文水动力过程，对于人类揭示自然奥秘、研究气候变化及适应策略、开发利用水及水能资源、发展航运、开展洪涝灾害管理、进行水生态环境保护等方面意义重大。

但随着计算机技术的发展，使用数值模型手段模拟水文水动力过程成为可能，相比于传统手段，数值模型可以在短时间构建完成并便捷计算不同条件下的结果，且可更加精细地描述物理过程，故应用范围越来越广泛。但传统的模拟在面对高分辨率的真实地形时其计算稳定性和精度仍存不足。因为高分辨率地形的复杂性导致地表水动力过程的复杂性，在计算中会导致异常水深和伪高流速等非物理现象，造成计算失稳和物质动量的不守恒。此外，采用高分辨率地形模拟，时间步长短和网格单元多会导致计算效率大为降低，对于管网排水效果的模拟，现有模拟方法跟多的采用其他等效替代方法，没有建立完整的水动力过程。

技术实现要素：

本发明的目的是提供一种基于显卡加速耦合管网的高精度城市雨洪模型，解决了现有技术中存在的传统雨洪模拟面对高分辨率的真实地形时其计算稳定性和精度不足的问题，实现高效高分辨率雨洪多过程耦合模拟。

本发明所采用的技术方案是，一种基于显卡加速耦合管网的高精度城市雨洪模拟方法，具体按照以下步骤实施：

步骤1、在计算机中读取计算区域的基础数据，包括地形数据dem、土地利用、曼宁系数、雨水井节点、管道数据、降雨数据和下渗数据；

步骤2、建立结构网格之间的拓扑关系以及网格与土地利用、曼宁的对应关系；

步骤3、依据雨水井节点坐标与地形网格坐标，依据就近原则建立雨水井节点与网格的对应关系；

步骤4、建立雨水井节点与管道的拓扑关系；

步骤5、将步骤1、2、3和4中的所有变量及对应关系在gpu中分配内存，并将上述变量的值复制到gpu变量中；

步骤6、采用godunov格式的有限体积法离散二维浅水方程，得到可以被计算机求解的降雨源项、通量项、底坡源项、摩阻源项和时间步进，然后在gpu上分别计算降雨源项、通量项、底坡源项、摩阻源项和时间步进；

步骤7、根据堰流公式和孔流公式计算地表与雨水井交换的流量；

步骤8、采用有限差分法求解一维圣维南方程组获得管道内非恒定水流入流流量和处理流量以及流速、水深等相关要素；

步骤9、雨水井节点水深、管道水深和对应地表网格水深更新计算；

步骤10、采用cudamemcpy函数将gpu上计算的结果值复制到主机内存中，并输出模拟的地表水深分布图、管网排水过程图。

步骤6中下渗量、通量项、底坡源项、摩阻源项和时间步进的计算方法分别对应green-ampt入渗模型、黎曼求解器、底坡通量法、半隐式法和二阶龙格-库塔法。

步骤6中，二维浅水方程具体为：

式中，t为时间，s；r为净雨速率，等于降雨强度减去入渗速率及管网排出水量；q为变量矢量包括水深h，两个方向上的单宽流量qx和qy；u、v为x、y方向上的流速；f、g为x、y方向上的通量矢量；s为源项矢量，包括降雨或下渗源项r、底坡源项及摩阻力源项；zb为河床底面高程，m；床面摩擦系数cf＝gn²/h^1/3，n为曼宁系数。

步骤6中，green-ampt入渗模型的方程具体为：

式中，fp为土壤下渗速率，cm/min；ks为饱和导水率，cm/min；d为土壤表层积水深度，cm；ls为湿润锋深度，cm；为毛管吸力，cm。

步骤7中，堰流公式和孔流公式计算过程为：

地表与雨水井的水量交换可分为两种情况，节点入、回流以及节点溢流。其交换量可采用下式计算：

式中，q为地表与雨水井交换的流量，m³/s；hj为节点水位，m；hw为该节点在地表对应网格处的水位水位，m；my为堰流系数，视情况取值，范围[0,1]；ap为雨水井节点面积，m²；b为雨水口宽度或节点周长，m；g为重力加速度，m/s²；z为地表高程，m；mk为孔口流量系数，视情况取值，范围[0,1]。

步骤8中采用扩散波法计算管道入流和出流流量，选用有限差分法求解的过程为：

扩散波方程如式(4)、式(5)所示。

式中，a为管道过水断面面积，m²；q为管道流量，m³/s；s为固定横截面沿流程的距离；t为时间，s；g为重力加速度；dz/ds为水力坡降；sf为摩阻比降；sf＝q|q|·n²/(a²·r^4/3)。

步骤9中节点水深更新计算方法为：

hⁿ＝h^n-1+4(q+qⁿin-qⁿout)dt/πd²(6)

式中，h^n-1为上一时间间隔内的节点水深，m；qⁿin和qⁿout分别为当前时间步长内的入流流量和出流流量，m³/s；dt为单位时间，s；q为节点的入流流量,m³/s，d为节点直径，m。

步骤9中管道水深更新计算方法为：

aⁿ＝(v^n-1+(qⁿin-qⁿout)dt)/l(7)

式中，qⁿin和qⁿout分别为当前时间步长内的管道的入流流量和出流流量，m³/s；v^n-1为上一时刻内管道内水的体积，m³；l为管长，m。

步骤9中雨水井对应的地表网格水深计算方法为：

zcouple-new＝zcouple-qdt/acell(8)

式中，zcouple-new为更新后的地表耦合网格水面高程，m；zcouple为之前地表耦合网格水面高程，m；acell为地表网格面积，m²。

本发明的特点在于，为解决传统模拟在计算稳定性、精度以及计算效率不足的问题，地表坡面汇流过程采用godunov格式的有限体积法求解二维浅水方程，使用黎曼求解器、底坡通量法、半隐式法和二阶龙格-库塔法计算通量项、底坡源项、摩阻源项和时间步长。

本方法下渗计算采用具有物理含义的green-ampt入渗模型，该入渗模型建立在毛细管吸水原理的基础上，在不考虑下渗动水头对下渗的影响，湿润锋之后的土体含水率均为饱和含水率，并且湿润锋的前后两侧土体含水率差异很大，水分的渗漏速度只有在上层毛细管孔隙达到饱和后才开始增加。

该方法同时耦合了地下管网模型，可采用全水动力方法模拟完整的降雨径流-管网排水-地表积涝过程，管网模型采用有限差分法求解扩散波，可准确模拟管网排水过程。

特点还在于集成了显卡(graphicalprocessingunit,gpu)并行加速技术，由于gpu拥有成百上千个计算核心，极大提高了模型计算效率，时效性强，可快速准确模拟较大范围高分辨率雨洪过程。

本发明的有益效果是，本发明的一种基于显卡加速耦合管网的高精度城市雨洪模拟方法，计算准确、用时较短，可实现较大范围的高分辨率雨洪过程模拟预测。

附图说明

图1是本发明一种基于显卡加速耦合管网的高精度城市雨洪模拟方法城市雨洪过程组成示意图；

图2是本发明一种基于显卡加速耦合管网的高精度城市雨洪模拟方法节点与网格的对应关系示意图；

图3是本发明一种基于显卡加速耦合管网的高精度城市雨洪模拟方法的gpu数据交换及求解过程流程图；

图4是本发明一种基于显卡加速耦合管网的高精度城市雨洪模拟方法的基于godunov格式的浅水方程有限体积数值方法求解过程图；

图5是本发明一种基于显卡加速耦合管网的高精度城市雨洪模拟方法green-ampt下渗模型示意图；

图6是本发明一种基于显卡加速耦合管网的高精度城市雨洪模拟方法节点入流流量求解示意图；

图7是本发明一种基于显卡加速耦合管网的高精度城市雨洪模拟方法节点溢流流量求解示意图；

图8是本发明一种基于显卡加速耦合管网的高精度城市雨洪模拟方法管道-节点计算示意图；

图9是本发明一种基于显卡加速耦合管网的高精度城市雨洪模拟方法节点水深求解示意图；

图10是本发明一种基于显卡加速耦合管网的高精度城市雨洪模拟方法管道水深求解示意图；

图11是本发明一种基于显卡加速耦合管网的高精度城市雨洪模拟方法本发明专利的某研究区域内涝峰值时刻积水分布结果图；

图12是本发明一种基于显卡加速耦合管网的高精度城市雨洪模拟方法本发明专利的某研究区域实测降雨排口流量过程及对比图。

具体实施方式

下面结合附图和实施例对本发明进行详细说明。

一种基于显卡加速耦合管网的高精度城市雨洪模拟方法，具体按照以下步骤实施：

步骤1，在计算机中读取计算区域的基础数据，主要包括地形数据(dem)、土地利用、曼宁系数、雨水井节点、管道数据、降雨数据和下渗数据等，城市雨洪过程组成示意图如图1所示；

步骤2、建立结构网格之间的拓扑关系以及网格与土地利用、曼宁对应关系；

步骤3、依据雨水井节点坐标与地形网格坐标，依据就近原则建立雨水井节点与网格的对应关系，网格与节点对应关系如图2所示；

步骤4中拓扑关系包括节点所连接的管道编号，管道前段节点编号和末端节点编号，以及节点和管道的属性。

步骤5中采用cuda库函数的cudamalloc进行变量分配，采用cudamemcpy函数进行变量值的复制，gpu数据交换及求解过程流程图如图3所示。

步骤6中基于godunov格式的浅水方程有限体积数值方法求解过程如图4所示。其中下渗量、通量项、底坡源项、摩阻源项和时间步进的计算方法分别对应green-ampt入渗模型、黎曼求解器、底坡通量法、半隐式法和二阶龙格-库塔法。

步骤6中二维浅水方程具体为：

式中，t为时间，s；r为净雨速率，等于降雨强度减去入渗速率及管网排出水量；q为变量矢量包括水深h，两个方向上的单宽流量qx和qy；u、v为x、y方向上的流速；f、g为x、y方向上的通量矢量；s为源项矢量，包括降雨或下渗源项r、底坡源项及摩阻力源项；zb为河床底面高程，m；床面摩擦系数cf＝gn²/h^1/3，n为曼宁系数。

步骤6中green-ampt入渗方程具体为：

式中，fp为土壤下渗速率，cm/min；ks为饱和导水率，cm/min；d为土壤表层积水深度，cm；ls为湿润锋深度，cm；为毛管吸力，cm。下渗模型示意图如图5所示。

步骤7中堰流公式和孔流公式计算过程为：

地表与雨水井的水量交换可分为两种情况，节点入、回流以及节点溢流，节点入流和节点溢流求解过程分别如图6和图7所示。其交换量可采用下式计算：

式中，q为地表与雨水井交换的流量，m³/s；hj为节点水位，m；hw为该节点在地表对应网格处的水位水位，m；my为堰流系数，视情况取值，范围[0,1]；ap为雨水井节点面积，m²；b为雨水口宽度或节点周长，m；g为重力加速度，m/s²；z为地表高程，m；mk为孔口流量系数，视情况取值，范围[0,1]。

步骤8中一维圣维南方程组描述的非恒定水流是一种浅水中的长波现象，称为动力波。忽略公式中的惯性项，所得到的波称为扩散波。扩散波可以较好地反映水流在管网中的运动状态，因此本专利采用扩散波法计算管道入流和出流流量，选用有限差分法求解。

扩散波方程如式(4)、式(5)所示。

式中，a为管道过水断面面积，m²；q为管道流量，m³/s；s为固定横截面沿流程的距离；t为时间，s；g为重力加速度；dz/ds为水力坡降；sf为摩阻比降；sf＝q|q|·n²/(a²·r^4/3)。管道-节点计算示意图如图8所示。

步骤9中节点水深更新计算方法为：

hⁿ＝h^n-1+4(q+qⁿin-qⁿout)dt/πd²(6)

式中，h^n-1为上一时间间隔内的节点水深，m；qⁿin和qⁿout分别为当前时间步长内的入流流量和出流流量，m³/s；dt为单位时间，s；q为节点的入流流量,m³/s，d为节点直径，m。节点水深求解示意图如图9所示。

步骤9中管道水深更新计算方法为：

aⁿ＝(v^n-1+(qⁿin-qⁿout)dt)/l(7)

式中，qⁿin和qⁿout分别为当前时间步长内的管道的入流流量和出流流量，m³/s；v^n-1为上一时刻内管道内水的体积，m³；l为管长，m。管道水深求解示意图如图10所示。

步骤9中雨水井对应的地表网格水深计算方法为：

zcouple-new＝zcouple-qdt/acell(8)

式中，zcouple-new为更新后的地表耦合网格水面高程，m；zcouple为之前地表耦合网格水面高程，m；acell为地表网格面积，m²。

步骤10、采用cudamemcpy函数将gpu上计算的结果值复制到主机内存中，并输出模拟的地表水深分布图、管网排水过程图。

图11为某研究区域内涝峰值时刻积水分布结果图，其中a、b、c、d、e和f共6处积水点信息与实测值基本一致，表明本模拟系统在城市内涝模拟方面的可靠性。

图12为某研究区域实测降雨排口流量过程及对比图，结果显示本系统管网排水过程与实测管网排水过程接近，并且相较与其他模型模拟结果更佳。

本发明的有益效果是，本发明的一种显卡加速的高精度城市雨洪模拟方法，计算准确、用时较短，可实现较大范围的高分辨率雨洪过程模拟预测。