考虑无机态物质的GPU加速水环境模拟方法

本发明属于河湖海域水环境模拟，具体涉及考虑无机态物质的gpu加速水环境模拟方法。

背景技术：

1、由于工业化和城镇化社会的飞速发展，产生了大量的工农业废水和居民生活污水。这些污水中含有丰富的氮、磷等物质，若处理不充分被直接或间接地排放到河湖海域以及库、渠水体中。长期如此，造成水体污染难以恢复，水生生物死亡，甚至对人类健康形成直接威胁。为有效开展水环境污染防治与治理工作，水环境模型受到很多学者的青睐。然而，早期简单的bod-do水质模型及修正的多组分水质模型均无法满足河流水环境的模拟需求，需要构建多维、多介质和动态模拟为特征的水环境模型以应对复杂的河流污染。但模型内各环境因子之间存在复杂的转化关系，导致计算复杂，模拟时间较长。基于此，构建一种考虑无机态物质的gpu加速水环境模拟方法，以各物质状态变量转化平衡方程为基本理论依据，精准、高效地模拟河湖海域水体内各环境变量的输移规律及其相互之间的转化关系，从而实现河湖海域复杂水体富营养化状况的精确模拟与有效评估。

技术实现思路

1、本发明的目的是提供考虑无机态物质的gpu加速水环境模拟方法，能够实现复杂河湖海域水体中各环境变量输移规律及其变量之间相互转化关系的精细化模拟，以精确、有效地评估河湖海域水体的富营养化状况，可为水环境污染的有效预测和防治提供一个重要的模拟工具。

2、本发明所采用的技术方案是，考虑无机态物质的gpu加速水环境模拟方法，具体按照以下步骤实施：

3、步骤1、利用arcgis处理地形散点数据形成高分辨率dem地形表征河湖海域的真实地形，作为地形数据的输入条件；

4、步骤2、通过矩形网格离散上述dem地形数据，并根据已知边界条件对相关参数进行赋值；

5、步骤3、构建水动力-物质输运耦合方程；

6、步骤4、将水动力模型与上述水动力-物质输运耦合方程进行耦合；

7、步骤5、将水动力-物质输运耦合方程组源项中生化反应项ssch，利用剖开算子法以各物质状态变量转化平衡方程进行逐一表述；

8、步骤6、采用改进静水重构法，以满足复杂网格地形上静水模拟中通量项与底坡源项的平衡性；

9、步骤7、利用muscl二阶精度和二阶龙格库塔法实现模型的时空二阶精度；

10、步骤8、输出各环境因子的输移过程分布图，以反映各环境因子在水体内的输移规律，并评估水体富营养化状况。

11、本发明的特点还在于，

12、步骤3中水动力-物质输运耦合方程表示如下：

13、

14、式中，cs为沿水深方向的平均物质浓度，单位为kg/m3；dx、dy分别为x、y方向的紊动扩散系数，单位为m2/s；源汇项中ssch可表示为生化反应项，单位为kg/m3；c0为沿水深方向的平均物质浓度，单位为kg/m3。

15、步骤4中二维水动力-物质输运耦合方程组的矢量形式表示如下：

16、

17、其中，

18、

19、式中，t为时间，单位为s；q为变量矢量，其中，h为水深，单位为m，qx和qy分别为x、y方向上的单宽流量，单位为m2/s；f、g分别为x、y方向上的通量矢量；u、v分别为x、y方向上的流速，单位为m/s；s为源项矢量，包括底坡源项、摩阻源项、温度扩散项和温度源项；zb为河床底面高程，单位为m；cf为床面摩擦系数，计算式为cf＝gn2/h1/3，其中n为曼宁系数；g为重力加速度，单位为m/s2；cs为沿水深平均的物质平均浓度，单位为kg/m3；dx、dy分别为x、y方向的紊动扩散系数，单位为m2/s；

20、源汇项中ssch表示为生化反应项，单位为kg/m3；c0为沿水深方向的平均物质浓度，单位为kg/m3。

21、步骤4具体按照以下步骤实施：

22、步骤4.1、控制单元网格i内，采用中心有限体积法对二维水动力-物质输运耦合控制方程组进行积分，基于高斯散度定理，将通量项面积积分转化为线积分，具体为：

23、

24、式中，ω为控制容积i的体积；г为控制容积i的边界；n为控制容积边界г所对应的外法线方向的单位向量；f(q)·n为控制容积边界г上的外法向对流项通量向量；s为源项向量，包括底坡源项、摩阻源项，扩散项和生化反应项；

25、通量f(q)·n线积分为：

26、

27、式中，ω为第i个控制单元网格的体积，单位为m3；г为第i个控制单元网格的边界；n为边界г所对应的外法线方向的单位向量；sb表示底坡源项通量；sf表示摩阻源项通量；k为控制单元网格边的编号；lk为第i个控制单元网格第k个边的边长；

28、对通量项采用hllc黎曼求解器进行求解：

29、关于波速计算如下：

30、

31、其中，中间变量h*和根据式计算，

32、

33、

34、f(q)·n的计算结果为：

35、

36、其中，控制体单元网格界面中间通量f*为，

37、

38、控制体单元网格界面通量中间波两侧，和为，

39、

40、式中，fl和fr分别为单元网格界面左、右侧的通量；f*、和分别为控制体单元网格界面通量中间波两侧的左波通量、中间波通量和右波通量；u||为沿交界面的切向速度，u||＝-uny+vnx，u||l和u||r分别为控制体单元网格交界面的左侧、右侧切向速度，单位为m3/s；nx和ny分别为沿x和y的指示方向；u和v分别为沿x和y的单宽流量，单位为m/s；变量q⊥＝[h,qxnx+qyny]；qx和qy分别为沿x和y的单宽流量，单位为m2/s；sl、sm、sr分别为左波速、中间波速、右波速，单位为m/s；u⊥为通过网格单元界面垂直的速度，单位为m/s；u⊥＝unx+vny；h*和为中间变量；g为重力加速度，取9.81，单位为m/s2；cs为沿水深平均的物质平均浓度，kg/m3；

41、步骤4.2、将扩散项作为源项进行处理，采用二阶中心差分法对扩散项进行离散求解，物质输运扩散项表示为：

42、

43、根据二阶中心差分法，对物质输运扩散项的一阶导数和二阶导数进行求解，对于任意二维单元网格中心(i,j)得到：

44、

45、从而，得到物质输运扩散项的离散形式为：

46、

47、步骤4.3、底坡源项采用底坡通量法求解，得到某一边f处的底坡通量的矢量形式为：

48、

49、式中，fsf(q)为某一边f处的底坡通量；nf为面f单位外法向量；nfx和nfy是nf与(面f单位外法向量)的组成部分，hl和zbl分别为网格单元中心处的水深和底部高程，与zml网格边上的水深和底部高程，g为重力加速度，单位为m/s2；

50、步骤4.4、摩阻源项采用显隐式摩阻处理法求解，计算公式为：

51、

52、其中，单宽流量更新为时间步长n+1时刻计算公式为：

53、

54、式中，sfx和sfy分别为sf在x和y方向的摩阻源项；cf为床面糙率系数，cf＝gn2/h1/3，n为曼宁系数，s/m1/3；h为水深，单位为m；qx、qy分别为x、y方向上平均流速，单位为m2/s；qxn+1、qyn+1分别为时间步长n+1时刻在x、y方向上的单宽流量，单位为m2/s；qxn、qyn分别为时间步长n时刻在x、y方向上的单宽流量，单位为m2/s；δt为时间步长，单位为s；δqx为时间步长n时刻和n+1时刻的x方向上单宽流量之差，即δqx＝qxn+1-qxn，单位为m2/s；为了避免分母为零，θ取1.0e-12。

55、步骤5具体按照以下步骤实施：

56、步骤5.1、溶解氧平衡包括溶解氧do、生化需氧量bod和化学需氧量cod；溶解氧do浓度包括水生植物呼吸、硝化、生物降解和底泥或泥沙耗氧量以及光合作用产氧量；

57、步骤5.2、氮元素循环过程包括氨氮nh3、硝酸盐氮no3和亚硝酸盐氮no2三种形态；

58、步骤5.3、磷循环中磷酸盐po4浓度包括水生植物和细菌微生物消耗量，有机物降解反应和异氧生物呼吸作用的释放量；

59、步骤5.4、以浮游植物的生长动力学为主线，以c:n:p三个主要元素的比例反映浮游植物即叶绿素α与环境中营养盐之间的转化关系。

60、步骤6具体按照以下步骤实施：

61、步骤6.1、基于静水重构法，干湿单元边界边f中心点m处的底坡高程取值为：

62、

63、式中，分别为干湿单元边界边f中心点m处左右两侧底坡高程，单位为m；

64、重构两侧水深，同时，使干湿界面处右侧水位与重构后底坡高程相同，即：

65、

66、式中，分别为干湿单元边界边f中心点m左、右侧水位，单位为m；；分别为重构后左、右侧水深，单位为m；zbm为干湿单元边界边f中心点m重构后底坡高程，单位为m；

67、若用水位表示，则重构后的水位为：

68、

69、式中，分别为重构后水位，单位为m；

70、步骤6.2、在上述基础上，重构界面两侧的流量：

71、

72、式中，分别为重构后在x、y方向左、右侧单宽流量，单位为m2/s；分别为干湿单元边界边f中心点m在x、y方向上左、右侧流速，单位为m/s；

73、步骤6.3、干湿边界单元边f中心点m处的底坡高程值应取该点重构后底部高程和左侧水位的最低值，即：

74、

75、步骤6.4、进一步修正um、vm，确保干湿单元边界边单调变化，修正公式如下：

76、

77、式中，qxm、qym分别为干湿单元边界边f中心点m重构后在x、y方向上的单宽流量，单位为m3/s；um、vm分别为在x、y方向的相应流速，单位为m/s；hm为重构后相应水深，单位为m；ε为允许公差，由ε＝t2得到，t为最大单元的面积；

78、在复杂网格上通过将干湿边界界面处水深变化值和水深共同作为辨别条件，使得满足复杂网格地形的二阶或更高阶精度格式的干湿边界求解问题，判别公式如下：

79、

80、式中：为干湿单元边界边f中心点m左侧水深，单位为m；为干湿单元边界边f中心点m左侧的坡高程，单位为m；hi为i单元网格水深，单位为m；为i单元网格左侧的底坡高程，单位为m；γ为网格单元边缘处水深重构值与单元形心处水深重构值之比，γ＜1，γ值越高，符合这一标准的网格单元就越多，εwd为水深容差，通常取εwd＝1×10-6。

81、步骤7具体按照以下步骤实施：

82、步骤7.1、采用muscl重构法计算单元网格内各变量的空间变化梯度；通过限制变量梯度获得限制变量梯度，基于限制梯度，插值出单元边界中点处的变量值，以实现muscl二阶重构。

83、步骤7.2、由两步runge-kutta方法获得二阶时间精度，在新时间步长中第i个控制单元网格的值为：

84、

85、中间变量值为

86、

87、式中，i为控制单元网格数；ω为第i个控制单元网格的体积，单位为m3；qn+1为时间步长n+1时刻的单宽流量，单位为m2/s，qn为时间步长n时刻的单宽流量，单位为m2/s；δt为时间步长，单位为s；s(qn)为时间步长n时刻的源项；整个求解过程中，引入gpu并行计算技术以提高计算效率；

88、步骤7.3、在空间步长δx下，时间步长δt的取值的计算公式如下：

89、

90、式中，δxmin表示网格单元中心到对应界面的最小距离，单位为m；u为流速，单位为m/s；g为重力加速度，取值9.81m/s2；h为水深，单位为m。

91、本发明的有益效果是，考虑无机态物质的gpu加速水环境模拟方法，针对现有水环境模型中由于各环境因子之间存在复杂的转化关系，导致的计算复杂，模拟时间较长等缺点，提出的考虑无机态物质的gpu加速水环境模拟方法，精准、高效地模拟河湖海域水体内各环境变量的输移规律及其相互之间的转化关系，从而实现河湖海域复杂水体富营养化状况的精确模拟与有效评估。