CN105893763A - 南方感潮河网区水流水质耦合模拟*** - Google Patents

Abstract

本发明提供一种南方感潮河网区水流水质耦合模拟***，分为数据输入、模型库和数据输出三大模块，模型库包括水流数学模型和水质数学模型。水流数学模型可以作为单独的模块使用，对输入的地形、取水工程数据与水文数据等进行分析计算，采用灵活的用户界面使得数据可以方便地进行用户需要的二次处理，并为水质数学模型提供了基础数据支撑。水质数学模型与水流数学模型有机结合，实现了在水流数学模型的基础上对大范围网河区河流水质的模拟，并能够对闸坝、取排水等进行灵活设置，加强了***的实用性及适用性。

Description

南方感潮河网区水流水质耦合模拟***

技术领域

本发明涉及水环境数学模型技术领域，更具体地，涉及一种南方感潮河网区水流水质耦合模拟***。

背景技术

随着城市化的快速发展和城市人口的高速增长，城市供水需求不断增加，加之工业废水、生活污水的大量排放，使得本就突显的水资源供需矛盾及水污染问题更加剧烈。如何科学合理的进行取水以及污染物排放成为水资源可持续利用与发展和水环境保护急需解决的问题。

水环境数学模型是将水动力问题和污染物在水体中的迁移转化用问题用数学方程进行描述，并在一定的定解条件下求解这些方程，从而模拟出实际问题的水流、水质状况的有效工具，可用于指导水环境工程实践并实现对水环境的规划管理和污染控制，在水环境问题的研究中发挥着越来越重要的作用。

目前，国内外应用范围较为广泛的主要水环境数学模型包括QUAL2K河流水质综合模型、MIKE水动力与水质模型、WASP水质模型、SMS综合性模型等。其中：1)QUAL2K是由美国环保局(USEPA)开发并不断修订和增强的河流综合水质模型。是一种一维稳态模型，适用于枝状河流模拟，沿程可有多个支流汇入和多个排污口，还可模拟河道水工建筑物对水质的影响。国内戴凌全、李华等以新安江大坝附近河段为例讨论了该模型主要参数在实际应用中的确定方法。2)MIKE是丹麦DHI公司开发的水动力学与水质模型，MIKE11用于一维水动力、水质、富营养化和泥沙输移计算及洪水预报等，适用于河流湖库及灌溉渠等。MIKE21用于二维水动力学、水质、富营养化、石油泄漏等计算，适用于河流、湖库、河口及海湾等。MIKE31用于三维水动力和水质模拟，适用于河流、湖库、河口、海洋等。国内有用该模型进行汉江、太湖等水体富营养化的研究。3)WASP是由美国环保局开发的水质模型程序，能够用于分析预测自然和人为污染的各种水质状况，可模拟水文动力学、河流一维不稳定流、湖泊和河口三维不稳定流、常规污染和有毒污染物的迁移转化。4)SMS是由美国BrighamYoung大学环境模型研究实验室开发的一个综合性模型。它包含一、二维有限单元模型、有限差分模型以及三维水动力学模型。可用来进行河流、河口、湖泊、海岸水动力和水质模拟。国内将该模型用于工程实践的有五里湖调水模拟、黄河口水流数值模拟等。

感潮河网不同于山区单向河流。感潮河网内部结构错综复杂、水流潮汐往复，造成方程组离散及其求解非常困难，这是多年来人们研究河网问题的一大难点，也因此导致许多成熟的水环境数学模型往往不再适用。且国内现有的相关技术主要对理论研究较多，对模型以及软件的开发缺少必要的投入，可用于实际操作的成果较少。另外，在模型的设计过程中忽略了复杂的河道地形、初边条件和潮水涨落等因素，导致模型的稳定性、收剑性较差且使用过于复杂；此外在少数的操作***中，主要针对单一河流或者小范围流域的水流水质模拟较为普遍，对于大范围复杂网河区水流水质模拟不适用或精度低；再一方面，国内现有开发***开发年代较早，没有考虑当下剧烈人类活动，如修筑水闸、大坝、桥梁、取水工程等对河流水动力过程的影响作用，较难满足当前新情况、新问题频出的水资源开发利用及水环境保护情势。

发明内容

本发明为克服上述现有技术所述的至少一种缺陷，提供一种南方感潮河网区水流水质耦合模拟***(SWRNQ)。其目的是通过水流数学模型模拟南方感潮河网区的水动力过程，并结合水质数学模型，对取水工程以及污染物排放等进行水量水质演算，为南方地区取用水和排污制定合理方案，达到水资源的可持续开发利用和管理，实现水环境保护目标。

为解决上述技术问题，本发明的技术方案如下：

南方感潮河网区水流水质耦合模拟***，所述***包括：

数据输入模块：用于输入数据，输入数据由水流模拟、水质模拟和水流水质模型耦合关系三部分所需数据文件组成，其中水流模拟所需文件包括河道流场边界文件，流场边界类型及流场河道节点文件，河道地形数据文件，取水位置和取水量文件；水质模拟所需文件包括浓度场河道首末节点文件，浓度场污染源文件，浓度场外河道地形走向判断因子文件；水流水质耦合关系所需文件为浓度场与流场的河道、节点对应关系文件；

模型库：包括水流数学模型和水质数学模型两大模块，且两者相互连接，有机结合；所述水流数学模型导入和存储河道地形数据、模型边界水文数据、取用水数据和模型河道糙率参数初始值，然后根据输入数据自动进行模型参数率定，最后使用率定后的模型进行工程取水论证，以期达到水资源的可持续开发和利用，并得到水质数学模型所需的基础数据；所述水质数学模型在水流水质模型耦合关系下通过水流数学模型提供的基础数据，结合输入的污染源位置及排污量数据、模型边界水质浓度数据和污染物降解系数等进行河道水质模拟演算，制定合理的排污方案；

数据输出模块：用于输出数据，数据输出文件由水流数学模型及水质模型两部分的输出结果组成，其中潮位和流量结果由水流数学模型输出，污染物浓度结果由水质数学模型输出。

在一种优选的方案中，所述水流数学模型以一维非恒定流圣维南方程组为基础，采用四点加权Preissmann隐式差分格式对方程进行离散，以提高模型稳定性及模拟精度；采用矩阵标识法对节点方程进行求解，以提高计算效率，实现节点任意编码，增强可扩充性、可移植性和通用性；根据南方感潮河网区河道水闸众多的现状可自行添加水闸，以更真实地反映河网区的水动力过程；对于河道糙率率定提供了一套适用于南方感潮河网区的调整规则。

在一种优选的方案中，所述水质数学模型采用适用于各种流动情况下求解河网对流输移问题隐式差分方程的一套计算公式及相应的河网计算方法，对于网河区某一河道，根据河道两端的水流方向的不同组合设置了4种流态，即：顺流向的流动、逆流向的流动、河道两端向中间的流动、河道中间任意河段向两端的流动；在河网计算中，采用了一种结合河道节点方程来补充水质方程的个数的求解方法。

本发明技术方案的有益效果是：

(1)***分为数据输入、模型库和数据输出三大模块，将地形、取水工程数据与水文数据等单独在水流数学模型进行分析处理，水质数学模型在水流数学模型分析处理的前提下进行水质模拟计算，这种前后继承关系的模块化设计实现了***的有机结合，确保了***的高效运行及稳定性。

(2)水流数学模型可以作为单独的模块使用，对输入的地形、取水工程数据与水文数据等进行分析计算，采用灵活的用户界面使得数据可以方便地进行用户需要的二次处理，并为水质数学模型提供了基础数据支撑。采用矩阵标识法对河网进行编码，适用性和可移植性大大提高。

(3)水质数学模型与水流数学模型有机结合，实现了在水流数学模型的基础上对大范围网河区河流水质的模拟，并能够对闸坝、取排水等进行灵活设置，加强了***的实用性及适用性。

附图说明

图1南方感潮河网区水流水质耦合模拟***结构图。

图2南方感潮河网区水流水质耦合模拟***运行示意图。

具体实施方式

下面结合附图和实施例对本发明的技术方案做进一步的说明。

实施例1

如图1-2所示，南方感潮河网区水流水质耦合模拟***，包括：

数据输入模块：用于输入数据，输入数据由水流模拟、水质模拟和水流水质模型耦合关系三部分所需数据文件组成，其中水流模拟所需文件包括河道流场边界文件(zqb.dat)，流场边界类型及流场河道节点文件(nrscb.dat)，河道地形数据文件(hdsj.dat)，取水位置和取水量文件(getwater.dat)；水质模拟所需文件包括浓度场河道首末节点文件(krc.dat)，浓度场污染源文件(scp.dat)，浓度场外河道地形走向判断因子文件(kib.dat)；水流水质耦合关系所需文件为浓度场与流场的河道、节点对应关系文件(knrcc.dat)；

模型库：包括水流数学模型和水质数学模型两大模块，且两者相互连接，有机结合；所述水流数学模型导入和存储河道地形数据、模型边界水文数据、取用水数据和模型河道糙率参数初始值，然后根据输入数据自动进行模型参数率定，最后使用率定后的模型进行工程取水论证，以期达到水资源的可持续开发和利用，并得到水质数学模型所需的基础数据；所述水质数学模型在水流水质模型耦合关系下通过水流数学模型提供的基础数据，结合输入的污染源位置及排污量数据、模型边界水质浓度数据和污染物降解系数等进行河道水质模拟演算，制定合理的排污方案；

数据输出模块：用于输出数据，数据输出文件由水流数学模型及水质模型两部分的输出结果组成，其中潮位和流量结果由水流数学模型输出，污染物浓度结果由水质数学模型输出。

下面对输入数据模块的各文件做具体说明：

(1)流场边界类型及流场河道节点文件

根据上下游流量水位边界不同，对不同边界进行识别，1为水位边界，0为流量边界，如表1所示，对河网区河道进行编码，先编内河道，再编外河道；对河道连结处进行节点设置、编码，先编内节点，再编外节点。按照流场河道的先后顺序，以及连接的首节点和末节点编写河道节点文件。如表2所示。

表1流场边界类型文件格式(nrcb.dat)

河道边界类型

0,0,0,0,0,0,0,0,0,1,1,1,1

1为水位边界，0为流量边界

表2流场河道节点文件格式(nrcb.dat)

首节点编号	末节点编号
		1
		N

(2)河道流场边界文件

输入外河道上边界的流量数据和下边界的水位数据，其中，边界数据可以是历时的流量潮位数据历时典型的潮位数据。如表3所示。

表3河道边界文件格式(zqb.dat)

(3)河道地形文件

对已编码的河道按顺序录入河道的断面数量、断面糙率、断面地形和断面之间的距离。如表4所示。

表4河道地形文件格式(hdsj.dat)

(4)取水口文件

根据已有取水口的信息，在模型中输入取水口的数量、位置、取水量等信息。如表5所示。

表5取水口文件格式(getwater.dat)

河道序号	断面序号	取水量(m3/s)

(5)浓度场与流场的河道、节点对应关系文件

根据浓度场河道的编号先后顺序对应的流场河道编号，以及根据浓度场节点的编号的先后对应的流场节点编号。如表6所示。

表6浓度场与流场的河道、节点对应关系文件文件格式(knrcc.dat)

(6)浓度场河道首末节点文件

按照浓度场河道的先后顺序，以及连接的浓度场中首节点和末节点编写浓度场河道节点文件。如表7所示。

表7浓度场河道首末节点文件(krc.dat)

(7)浓度场污染源文件

根据已有排污口的信息，在模型中按照排污口的顺序输入位置、氨氮排放量、COD排放量等信息。如表8所示。

表8浓度场污染(scp.dat)

河道序号	断面序号	氨氮排放量(kg/year)	COD排放量(kg/year)

(8)浓度场外河道地形走向判断因子

对于浓度场外河道断面编号一致是从外向里编号，这与流场中河道断面编号有可能不一致，1为走向不一致，2为走向一致。

表9流场边界类型文件格式(kib.dat)

下面对模型库做具体说明：

***中水流数学模型采用圣维南方程组为描述河道的水流运动的基本方程，采用三级联合解法。采用四点隐式差分格式离散微段控制方程组，求解时，将微段方程自相消元，得到河道首、末断面流量由首、末断面水位表示的河段方程，再利用节点水量平衡方程，得到以节点水位为基本未知量的节点方程组；针对节点方程组系数矩阵的高稀疏及主对角线占优的特点，采用矩阵标识法求解，解出节点水位后，回代到各河道，可求出各计算微段的流量和水位。

1)控制方程

描述河道非恒定流的基本方程为一维圣维南方程组：

连续方程

动量方程

式中，Z为断面水位；Q为流量；A为河道过水面积；g为重力加速度；B为过水宽度；q为旁侧入流流量；R为水力半径；c为谢才(Chezy)系数；x、t是位置和时间坐标。

2)差分格式

采用四点加权Preissmann隐式差分格式圣维南方程组式(1)和式(2)，若以F代表流量Q和水位Z，针对图中M点建立差分因子，则任意函数F及其偏导数的离散形式为：

$\left\{\begin{array}{c}F=\frac{1}{2}(F_{j+1}^n+F_j^n)\\ \frac{\∂F}{\∂x}=\mathrm{\θ}\frac{F_{j+1}^{n+1}-F_j^{n+1}}{\mathrm{\Δ}x}+(1-\mathrm{\θ})\left(\frac{F_{j+1}^n-F_j^n}{\mathrm{\Δ}x}\right)\\ \frac{\∂F}{\∂t}=\frac{F_{j+1}^{n+1}+F_j^{n+1}-F_{j+1}^n-F_j^n}{2\mathrm{\Δ}t}\end{array}\right.---\left(3\right)$

式中，θ为加权系数，一般取0.5-1.0。

3)微段方程

3.1)差分方程

设河道共有m个断面，则有m-1个微段，首断面编号为1，磨端面编号为m。按照(3)式所描述的简化四点线性隐式差分，代入连续方程(1)可得：

$\frac{B_{j+\frac{1}{2}}^n}{2\mathrm{\Δ}t}(Z_{j+1}^{n+1}-Z_{j+1}^n+Z_j^{n+1}-Z_j^n)+\mathrm{\θ}\left(\frac{Q_{j+1}^{n+1}-Q_j^{n+1}}{{\mathrm{\Δx}}_j}\right)+(1-\mathrm{\θ})\left(\frac{Q_{j+1}^n-Q_j^n}{{\mathrm{\Δx}}_j}\right)=q_{j+\frac{1}{2}}---\left(4\right)$

假设时段初的变量均为已知，整理得：

$-Q_j^{n+1}+Q_{j+1}^{n+1}+c_j{Z}_j^{n+1}+c_j{Z}_{j+1}^{n+1}=D_j---\left(5\right)$

其中：

$\left\{\begin{array}{c}{c}_j=\frac{B_{j+\frac{1}{2}}^n{\mathrm{\Δx}}_j}{2\mathrm{\Δ}t\mathrm{\θ}}\\ D_j=\frac{q_{j+\frac{1}{2}}{\mathrm{\Δx}}_j}{\mathrm{\θ}}-\frac{1-\mathrm{\θ}}{\mathrm{\θ}}(Q_{j+1}^n-Q_j^n)+c_y(Z_{j+1}^n+Z_j^n)\\ B_{j+\frac{1}{2}}^n=\frac{B_j^n+B_{j+1}^n}{2}\\ q_{j+\frac{1}{2}}=\frac{q_j^n+q_{j+1}^n}{2}\end{array}\right.---\left(6\right)$

对于动力方程(2)则可整理为

其中：

由曼宁公式则

为书写方便，忽略上标n+1，可把式的任一微段差分方程写为：

$\left\{\begin{array}{c}-Q_j+Q_{j+1}+c_j{Z}_j+c_j{Z}_{j+1}=D_j\\ E_j{Q}_j+G_j{Q}_{j+1}-F_j{Z}_j+F_j{Z}_{j+1}={\mathrm{\Φ}}_j\end{array}\right.---\left(9\right)$

其中c_j、D_j、E_j、F_j、G_j、φ_j均由初值计算，所以方程组为常系数线性方程组。对一条有m-1个微段的河道，有2(m-1+1)个未知量，可以列出2(m-1)个方程，加上河道两端的边界条件，形成封闭的代数方程组：

3.2)追赶方程

外河道(单一河道)的追赶方程

上边界为流量边界条件，追赶方程为：

$\begin{array}{c}\begin{array}{ccc}{Q}_1=P_1-V_1{Z}_1& P_1=Q_1\left(t\right),& V_1=0\end{array}\\ \{\begin{array}{c}{Z}_j=S_{j+1}-T_{j+1}{Z}_{j+1}\\ Q_{j+1}=P_{j+1}-V_{j+1}{Z}_{j+1}\end{array},j=1,2,\mathrm{......},m-1\end{array}---\left(11\right)$

其中：

$\left\{\begin{array}{ccc}{S}_{j+1}=\frac{G_j{y}_3-y_4}{y_1{G}_j+y_2},& T_{j+1}=\frac{G_j{C}_j-F_j}{y_1{G}_j+y_2}& \\ P_{j+1}=y_3-y_1{S}_{j+1},& V_{j+1}=C_j-y_1{T}_{j+1}& \\ y_1=V_j+C_j,& y_2=F_j+E_j{V}_j,& y_3=D_j+P_j,y_4={\mathrm{\Φ}}_j-E_j{P}_j\end{array}\right.---\left(12\right)$

下边界为水位，即Zm已知，末断面流量Q_m＝P_m-V_mZ_m，用Zm回代到方程组(11)，按j＝m-1,m-2,……,1的顺序求出微段断面的流量Qj及水位Zj，当j＝1时，求出Q2、Z1，而Q1为已知的上边界。

内河道的追赶方程

对于内河道，首、末断面的边界条件均为未知，由于没有端点边界条件可供利用，单一河道递推的方法不能适用。

内河道计算采用的思路为：以内节点水位为基本未知量，利用追赶方程，用逆推法(j＝m-1,m-2,……,1)、顺推法(j＝2,3,……,m)得到各河道各断面的流量用首、末断面水位表示的两个表达式，那么首、末断面的流量也可表示为首、末断面水位的唯一表达式。根据stoker条件，首、末断面的水位即与之相连节点的水位，再由水量平衡公式，得到节点水位平衡方程，由此得到各节点水位方程构成的方程组，从而求得内节点的水位，再求得河道首、末断面流量。回代到各微段方程，可求得各计算断面水位与流量。

逆推法：以方程组(10)为基础方程，从第m-1个微段方程起，用逆推法相互消元，把内河道中任一断面的流量表达成本断面水位和末断面水位的线性函数：

Q_j＝α_j+βZ_j-ξ_jZ_m，(j＝m-1,m-2,……,2,1) (13)

对于上式(13)，当j＝m-1时，有单独递推系数的起始值：

$\left\{\begin{array}{c}{\mathrm{\α}}_{m-1}=\frac{{\mathrm{\Φ}}_{m-1}-G_{m-1}{D}_{m-1}}{G_{m-1}+E_{m-1}}\\ {\mathrm{\β}}_{m-1}=\frac{C_{m-1}{G}_{m-1}+F_{m-1}}{G_{m-1}+E_{m-1}}\\ {\mathrm{\ξ}}_{m-1}=\frac{C_{m-1}{G}_{m-1}-F_{m-1}}{G_{m-1}+E_{m-1}}\end{array}\right.---\left(14\right)$

其余的递推公式为：

$\left\{\begin{array}{c}{\mathrm{\α}}_j=\frac{y_1({\mathrm{\Φ}}_j-{\mathrm{\α}}_{j+1}{G}_j)-y_2(D_j-{\mathrm{\α}}_{j+1})}{y_1{E}_j+y_2}\\ {\mathrm{\β}}_j=\frac{y_2{C}_j+y_1{F}_j}{y_1{E}_j+y_2}\\ {\mathrm{\ξ}}_{m-1}=\frac{{\mathrm{\ξ}}_{j+1}(y_2-y_1{G}_j)}{y_1{E}_j+y_2}\\ y_1=C_j+{\mathrm{\β}}_{j+1},y_2=G_j{\mathrm{\β}}_{j+1}{F}_j\end{array}\right.,(j=m-2,m-3,\mathrm{......},1)---\left(15\right)$

顺推法：以方程组(10)为基本方程，从第1河段微段起，用顺推法进行相互消元，可把河道任一断面的流量表达成本断面水位与首段面的水位线性函数：

Q_j＝θ_j+ηZ_j-γ_jZ₁，(j＝2,3,……,m) (16)

对于上式(16)，当j＝2时，有单独递推系数的起始值：

$\left\{\begin{array}{c}{\mathrm{\θ}}_2=\frac{E_1{D}_1+{\mathrm{\Φ}}_1}{E_1+G_1}\\ {\mathrm{\η}}_2=-\frac{C_1{E}_1+F_1}{E_1+G_1}\\ {\mathrm{\γ}}_2=\frac{F_1-C_1{E}_1}{E_1+G_1}\end{array}\right.---\left(17\right)$

其他系数的递推公式为：

$\left\{\begin{array}{c}{\mathrm{\θ}}_j=\frac{y_2(D_{j-1}+{\mathrm{\θ}}_{j-1})-y_1({\mathrm{\Φ}}_{j-1}-E_{j-1}{\mathrm{\θ}}_{j-1})}{y_2-G_{j-1}{y}_1}\\ {\mathrm{\η}}_j=\frac{y_1{D}_{j-1}-y_2{C}_{j-1}}{y_2-G_{j-1}{y}_1}\\ {\mathrm{\γ}}_j=\frac{{\mathrm{\γ}}_{j-1}(y_2-y_1{E}_{j-1})}{y_2-G_{j-1}{y}_1}\\ y_1=C_{j-1}-{\mathrm{\η}}_{j-1},y_2=E_{j-1}{\mathrm{\η}}_{j-1}-F_{j-1}\end{array}\right.,(j=3,4,\mathrm{......},m)---\left(18\right)$

当首、末断面水位求得后，联解(13)、(16)，对于同一断面上的流量有：

$\left\{\begin{array}{c}{Q}_j={\mathrm{\α}}_j+{\mathrm{\β}}_j{Z}_j+{\mathrm{\ζ}}_j{Z}_m\\ Q_j={\mathrm{\θ}}_j+{\mathrm{\η}}_j{Z}_j+{\mathrm{\γ}}_j{Z}_1\end{array}\right.,(j=2,3,\mathrm{......},m-1)---\left(19\right)$

可得：

$Z_j=\frac{{\mathrm{\θ}}_j-{\mathrm{\α}}_j+{\mathrm{\γ}}_j{Z}_1-{\mathrm{\ζ}}_j{Z}_m}{{\mathrm{\β}}_j-{\mathrm{\η}}_j}---\left(20\right)$

求出Zj后，代入(19)任一方程，即可求出Qj。

4)河段方程

对于外河道，末断面流量为：

Q_m＝f_m(Z_m) (21)

对于内河道：当j＝1时，由式(13)得首断面流量：

Q₁＝α₁+β₁Z₁+ξ₁Z_m (22)

当j＝m时，由方程式(15)得：

Q_m＝θ_m+η_mZ_m+γ_mZ₁ (23)

(21)、(22)、(23)式即形成河段首、末断面以水位和流量为状态变量的河段方程，河网每一河道都能得到上面的河段方程，首、末断面的流量分别用首、末断面的水位Z₁、Z_m表示，外河道末断面流量只用末断面水位Z_m表示，而Z₁、Z_m即首、末节点的水位。

5)汊点方程

汊点是河段的交汇点，在汊点上，水流必须满足水流连续(即水量守恒)及动量守恒条件。利用汊点相容方程和边界方程，消去河段首、末断面的某个状态变量(流量或水位)，形成节点水位或流量的汊点方程组。

由于河道数比汊点数要多，一个河道有首、末断面，如利用水位相等(简化动量守恒条件)消去水位变量，形成流量节点方程组，则方程组阶数较大，求解不便。而利用水量守恒原理，消去流量变量，一个节点可形成与该节点相邻的节点水位为未知量的线性代数方程，形成水位节点方程则阶数较小，求解较方便。

(1)水流连续条件

设出流为正，入流为负，(也可设入流为正，出流为负)，进出每一节点的流量必须与该节点内实际水量的增减率相平衡，可表示为：

$\underset{i=1}{\overset{l}{\Σ}}{Q}_i^k=\frac{\∂{\mathrm{\Ω}}_k}{\∂t}=A_k\frac{Z_k^{n+1}-Z_k^n}{\mathrm{\Δ}t}---\left(24\right)$

式中，l为汊点相连河段数，k为汊点号，Ω_k表示第k汊点的蓄水量，A_k为调蓄节点的蓄水面积， 分别为调蓄节点时段末、时段初的水位。

若调蓄节点面积很小，则可忽略水位变化引起调蓄节点水体积变化，即该汊点概化为一几何点，式(24)可简化为：

$\underset{i=1}{\overset{l}{\Σ}}{Q}_i^k=0---\left(25\right)$

(2)能量衔接条件

如果汊点可概化为一个几何点，出入该汊点的水流平缓，不存在水位突变的情况，则各节点相连汊道的水位应相等，等于该点的平均水位，即：

如果各断面的过水面积相差悬殊，流速有较明显的差别，当略去汊点的局部损耗时，得伯努利(Bernouli)方程，各断面之间的能量水头应相等，即：

$Z_1+\frac{{\mathrm{\μ}}_1^2}{2g}=Z_2+\frac{{\mathrm{\μ}}_2^2}{2g}---\left(27\right)$

也可以写为：

aZ₁+bQ₁+cZ₂+bQ₂＝0 (28)

当汊点可概化为一几何点，由式(25)、(26)构成的条件称为Stoker条件，这是Stoker1957年提出的方法，目前在许多实际应用中都采用这种简化的方法。

在假设(26)成立前提下，把每一河道首、末断面与节点水位关系式(21)、(22)、(23)代入(25)，得到与节点i相邻的节点水位为未知变量的线性代数方程：

f_i(Z_i,j)＝0 (29)

其中Z_i,j为节点i相邻节点水位的集合。

对河网每一个节点，都建立上述的节点水位方程，形成以河网节点水位为基本未知变量的线性方程组：

AZ＝R (30)

其中：为系数矩阵，Z＝[Z₁ Z₂ … Z_n]^T为节点水位列阵，R＝[r₁ r₂ … r_n]^T为右端项列阵。n为节点总数。

当i,j不是相邻节点时，或节点i与节点j相连，但连接节点i的流量方向并不是流向节点i，则a_ij＝0，AZ＝R为带型线性方程组。

采用矩阵标识法对节点方程组进行求解，其求解的基本思想是：根据节点水位方程系数矩阵的高稀疏性，对矩阵非零元素进行代码标识。按照代码指示，把非零元素用一维数组存贮，排除零元素，节约内存。求解时，由代码指示，只对非零元素进行运算，从而大大提高方程组求解计算的效率。

6)模型建立及求解

a.数据准备

节点编码：对内节点进行统一编码，编码顺序可以任意。

河道编码：分别对内河道与外河道统一编码，编码顺序可以任意。

断面编码：对于内河道，以首节点对应首断面，末节点对应末断面，断面编号从首断面向末断面递增，代表河道的计算流向；对于外河道，以外节点向内节点递增为原则。

外河道边界信息：按实际的边界条件类型确定(有水位边界，流量边界或水位流量关系边界)，若为水位边界，则边界条件资料为水位过程线，若为流量边界，则边界条件资料为流量过程线。

河道的计算信息：根据节点编码和断面编码，确定内河道的首节点和末节点号，确定外河道的末节点号，确定河道首断面号和末断面号。

b.计算步骤

给定初始条件：如各断面的初始水位、流量等

时层循环开始

①计算外河道的边界条件，系数矩阵初始化；

②对可调蓄节点，将蓄水量的变化表达成流量与水位的线性关系，叠加到相应节点水位方程；

③按外河道的顺序，依次计算外河道的追赶系数，将末断面的流量水位关系叠加到末节点水位方程；

④按内河道的顺序计算内河道的追赶系数，将首末断面的流量方程分别叠加到首末节点水位方程；

⑤用矩阵标识法求解节点水位方程；

⑥由节点水位回代到河段、微段方程，求出各计算断面的水位和流量。

转下一时层循环。

水质数学模型

南方河网地区，河道纵横交错，水流流向顺逆不定，污染物也随之沿河向不同方向输移，水流中输送的物质浓度属非稳态非均匀变化等特点。本***采用一套适用于各种流动情况下求解河网对流输移问题隐式差分方程的一套计算公式，及相应的河网计算方法。

1)水质基本方程

采用一维对流扩散输移方程：

河道方程：

节点方程：

若忽略节点的调蓄面积，即Ω＝0，则式(31)可写为：

$\underset{l=1}{\overset{NL}{\Σ}}{\left(QC\right)}_{l,j}=0_j---\left(32\right)$

式中：Q，Z是流量及水位；A是河道断面积；E_X是纵向分散系数；C是水流输送的物质浓度；Ω是河道节点的水面面积；j是节点编号；i是与节点j相联接的河道编号；S_C是与输送物质浓度有关的衰减项，对COD和氨氮可写为S_C＝K_dAC，K_d是衰减因子；S是外部的源或汇项。

2)水质模型求解方法

对基本方程(30)，用隐式差分迎风格式将微分方程离散。以顺流向情况的差分为例，式中的时间项采用前差分，对流项采用迎风差分，扩散项采用中心差分格式，得：

$\left\{\begin{array}{c}\frac{\∂\left(AC\right)}{\∂t}=\frac{{\left(AC\right)}_i^{k+1}-{\left(AC\right)}_i^k}{\mathrm{\Δ}t},\\ \frac{\∂\left(QC\right)}{\∂x}=\frac{{\left(QC\right)}_i^{k+1}-{\left(QC\right)}_{i-1}^{k+1}}{{\mathrm{\Δx}}_{i-1}},\\ \frac{\∂}{\∂x}\left({\mathrm{AE}}_x\frac{\∂C}{\∂x}\right)=\[\frac{{\left({\mathrm{AE}}_x\right)}_i^{k+1}{C}_{i+1}^{k+1}-{\left({\mathrm{AE}}_x\right)}_i^{k+1}{C}_i^{k+1}}{{\mathrm{\Δx}}_i}\\ -\frac{{\left({\mathrm{AE}}_x\right)}_{i-1}^{k+1}{C}_i^{k+1}-{\left({\mathrm{AE}}_x\right)}_{i-1}^{k+1}{C}_{i-1}^{k+1}}{{\mathrm{\Δx}}_{i-1}}\]\frac{1}{{\mathrm{\Δx}}_{i-1}},\\ S_c-S={\stackrel{\‾}{K}}_{d,i-1}^{k+1}{\left(AC\right)}_i^{k+1}-{\stackrel{\‾}{S}}_{i-1}^{k+1},\end{array}\right.---\left(33\right)$

对于逆流向情况可得到类似的结果，式中表示河段平均值，上角标k是时段的初值，k+1是时段末值，为书写方便，下文中凡出现时段末值，都省略写上标。

考虑到河网中流向顺逆不定，离散基本方程时，需要引入流向调节因子r_c及r_d，将顺、逆流向的离散方程统一到同一方程中，经整理后得

a_iC_i-1+b_iC_i+c_iC_i+1＝Z_i (i＝1,2,…,n) (34)

式中：a_i,b_i,c_i是系数；C_i是i断面时段末的浓度；n是某河道的断面数。对于一般断面(i＝2,Λ,n-1)有：

$\left\{\begin{array}{c}{a}_i=-\[r_{c1}{D}_{11}+r_{d1}{D}_{21}+F_{c1}\]\mathrm{\Δ}t/V,\\ b_i=\[r_{c1}{D}_{11}+r_{c2}{D}_{22}+r_{d1}{D}_{21}+r_{d2}{D}_{32}+F_{c2}-F_{d2}\]\mathrm{\Δ}t/V\\ +\[r_{c1}{\stackrel{\‾}{K}}_{d,i-1}+r_{d2}{\stackrel{\‾}{K}}_{d,i}\]\mathrm{\Δ}t+1.0,\\ c_i=-\[r_{c2}{D}_{22}+r_{d2}{D}_{32}-F_{d3}\]\mathrm{\Δ}t/V\\ z_i={\mathrm{\α}}_i{C}_i^k+\[r_{c1}{\stackrel{\‾}{S}}_{i-1}{\mathrm{\Δx}}_{i-1}+r_{d2}{\stackrel{\‾}{S}}_i{\mathrm{\Δx}}_i\]\mathrm{\Δ}t/V,\end{array}\right.---\left(35\right)$

对于首断面(i＝1)有：

$\left\{\begin{array}{c}{a}_1=0,\\ b_1=\[r_{d2}{D}_{32}-F_{d2}\]\mathrm{\Δ}t/V_2+r_{d2}{\stackrel{\‾}{K}}_{d,1}\mathrm{\Δ}t-1.0,\\ c_1=-\[r_{d2}{D}_{32}-F_{d3}\]\mathrm{\Δ}t/V_2,\\ z_1={\mathrm{\α}}_1{C}_1^k+r_{d2}{\stackrel{\‾}{S}}_1{\mathrm{\Δx}}_1\mathrm{\Δ}t/V_2,\end{array}\right.---\left(36\right)$

对于末断面(i＝n)有：

$\left\{\begin{array}{c}{a}_n=-\[r_{c1}{D}_{11}+F_{c1}\]\mathrm{\Δ}t/V_1,\\ b_n=\[r_{c1}{D}_{11}+F_{c2}\]\mathrm{\Δ}t/V_1+r_{c1}{\stackrel{\‾}{K}}_{d,n-1}\mathrm{\Δ}t+1.0,\\ c_n=0,\\ z_n={\mathrm{\α}}_n{C}_n^k+r_{c1}{\stackrel{\‾}{S}}_{n-1}{\mathrm{\Δx}}_{n-1}\mathrm{\Δ}t/V_1,\end{array}\right.---\left(37\right)$

其中：

$\left\{\begin{array}{cc}{V}_1={\mathrm{\Δx}}_{i-1}(A_{i-1}+A_i)/2,& V_2={\mathrm{\Δx}}_i(A_i+A_{i+1})/2,\\ V=r_{c1}{V}_1+r_{d2}{V}_2,& {\mathrm{\α}}_i=\frac{A_i^k}{A_i}\end{array}\right.---\left(38\right)$

$\left\{\begin{array}{cc}{D}_{11}={\left({\mathrm{AE}}_x\right)}_{i-1}/{\mathrm{\Δx}}_{i-1},& D_{22}={\left({\mathrm{AE}}_x\right)}_i/{\mathrm{\Δx}}_i,\\ D_{21}={\left({\mathrm{AE}}_x\right)}_i/{\mathrm{\Δx}}_{i-1},& D_{32}={\left({\mathrm{AE}}_x\right)}_{i+1}/{\mathrm{\Δx}}_i,\end{array}\right.---\left(39\right)$

$\left\{\begin{array}{cc}{F}_{C1}=(Q_{i-1}+{\mathrm{Qa}}_{i-1})/2,& F_{c2}=(Q_i+{\mathrm{Qa}}_i)/2,\\ F_{d2}=(Q_i-{\mathrm{Qa}}_i)/2,& F_{d3}=(Q_{i+1}-{\mathrm{Qa}}_{i+1})/2,\end{array}\right.---\left(40\right)$

上两式中的各个变量Qa是相应于流量Q的绝对值。

式(34)即是由n个方程组成的线性隐式差分方程组。若差分方程组的求解分单一河道的求解和节点方程的求解。

1、单一河道的求解

对于网河区某一河道，根据河道两端的水流方向的不同组合存在4种流态，即：顺流向的流动、逆流向的流动、河道两端向中间的流动、河道中间任意河段向两端的流动。随着流向的顺逆变化，在不同的时间和空间，方程组(34)中的各项系数将自动按式(35)～(37)取各种不同的数值，而且方程组的组合也不相同，因此需要针对不同的流动情况采用不同的方法进行求解。

顺流向的流动

由首断面至末断面的流动，可得到

$\left\{\begin{array}{c}{C}_1=C_1\left(t\right),\\ C_2={\mathrm{Wn}}_2+{\mathrm{On}}_2{C}_1+{\mathrm{Gn}}_2{C}_3,\\ C_3={\mathrm{Wn}}_3+{\mathrm{On}}_3{C}_1+{\mathrm{Gn}}_3{C}_4,\\ \mathrm{......}\\ C_i={\mathrm{Wn}}_i+{\mathrm{On}}_i{C}_1+{\mathrm{Gn}}_i{C}_{i+1},\\ \mathrm{......}\\ C_{n-1}={\mathrm{Wn}}_{n-1}+{\mathrm{On}}_{n-1}{C}_1+{\mathrm{Gn}}_{n-1}{C}_n,\\ C_n={\mathrm{Wn}}_n+{\mathrm{On}}_n{C}_1.\end{array}\right.---\left(42\right)$

式中的系数由下列各式所得

$\left\{\begin{array}{c}{\mathrm{Wn}}_1=0,{\mathrm{On}}_1=1,{\mathrm{Gn}}_1=0,\\ {\mathrm{Wn}}_i=\frac{z_i-a_i{\mathrm{Wn}}_{i-1}}{b_i+a_i{\mathrm{Gn}}_{i-1}},\\ {\mathrm{On}}_i=-\frac{a_i{\mathrm{On}}_{i-1}}{b_i+a_i{\mathrm{Gn}}_{i-1}},(i=2,\mathrm{...},n-1)\\ {\mathrm{Gn}}_i=-\frac{c_i}{b_i+a_i{\mathrm{Gn}}_{i-1}},\\ {\mathrm{Wn}}_n=\frac{z_n-a_n{\mathrm{Wn}}_{n-1}}{b_n+a_n{\mathrm{Gn}}_{n-1}},\\ {\mathrm{On}}_n=-\frac{a_n{\mathrm{On}}_{n-1}}{b_n+a_n{\mathrm{Gn}}_{n-1}},\\ {\mathrm{Gn}}_n=0.\end{array}\right.---\left(43\right)$

当首断面端的入流过程Q₁(t)及C₁(t)已知时，再根据各断面的流量、浓度初始条件，利用上面的公式可求得各断面的浓度。

逆流向的流动

由末断面至首断面的流动，可得到

$\left\{\begin{array}{c}{C}_1={\mathrm{Wl}}_1+{\mathrm{Ol}}_1{C}_n,\\ C_2={\mathrm{Wl}}_2+{\mathrm{Ol}}_2{C}_n+{\mathrm{Gl}}_2{C}_1,\\ \mathrm{......}\\ C_i={\mathrm{Wl}}_i+{\mathrm{Ol}}_i{C}_n+{\mathrm{Gl}}_i{C}_{i-1},\\ \mathrm{......}\\ C_{n-1}={\mathrm{Wl}}_{n-1}+{\mathrm{Ol}}_{n-1}{C}_n+{\mathrm{Gl}}_{n-1}{C}_{n-2},\\ C_n=C_n\left(t\right).\end{array}\right.---\left(44\right)$

其中：

$\left\{\begin{array}{c}{\mathrm{Wl}}_n=0,{\mathrm{Ol}}_n=1,{\mathrm{Gl}}_n=0,\\ Wli=\frac{z_i-c_i{\mathrm{Wl}}_{i+1}}{b_i+c_i{\mathrm{Gl}}_{i+1}},\\ {\mathrm{Ol}}_i=-\frac{c_i{\mathrm{Ol}}_{i+1}}{b_i+c_i{\mathrm{Gl}}_{i+1}},(i=n-1,\mathrm{...},2)\\ {\mathrm{Gl}}_i=-\frac{a_i}{b_i+c_i{\mathrm{Gl}}_{i+1}},\\ {\mathrm{Wl}}_1=\frac{z_1-c_1{\mathrm{Wl}}_2}{b_i+c_i{\mathrm{Gl}}_{i+1}},\\ {\mathrm{Ol}}_1=-\frac{c_1{\mathrm{Ol}}_2}{b_1+c_1{\mathrm{Gl}}_2},\\ {\mathrm{Gl}}_1=0.\end{array}\right.----\left(45\right)$

当末断面端的入流过程Q_n(t)及C_n(t)已知时，再根据各断面的流量、浓度初始条件，利用上面的公式可求得各断面的浓度。

河道两端向中间的流动

首、末断面都是入流的流动，可得到

当河道两端的入流过程Q₁(t)、C₁(t)及Q_n(t)、C_n(t)已知时，再根据初始条件，并利用式(42)或式(44)计算各系数，再用式(44)、式(45)即可进行求解。

河道中间任意河段向两端的流动

用通常的追赶法进行求解。

追的过程即把某一断面(i＝1,2,Λ,n)的浓度值表示为下一断面浓度值的函数：

$\left\{\begin{array}{c}{C}_i={\mathrm{Wd}}_i+{\mathrm{Gd}}_i{C}_{i+1}={\mathrm{Wh}}_i,(i=1,\mathrm{\Λ},n-1)\\ C_n={\mathrm{Wd}}_n.\end{array}\right.---\left(47\right)$

其中：

$\left\{\begin{array}{c}{\mathrm{Wd}}_1=z_1/b_1,\\ {\mathrm{Gd}}_1=-c_1/b_1,\\ {\mathrm{Wd}}_i=(z_i-a_i{\mathrm{Wd}}_{i-1})/(b_i+a_i{\mathrm{Gd}}_{i-1}),(i=2,\mathrm{...},n)\\ {\mathrm{Gd}}_i=-c_i/(b_i+a_i{\mathrm{Gd}}_{i-1}).\end{array}\right.---\left(48\right)$

赶的过程是从末断面逆序推倒其他各断面的浓度值，即：

$\left\{\begin{array}{c}{C}_n={\mathrm{Wd}}_n\\ C_i={\mathrm{Wd}}_i+{\mathrm{Gd}}_i{C}_{i+1}={\mathrm{Wh}}_i(i=n-1,n-2,\mathrm{...},1)\end{array}\right.---\left(49\right)$

将以上四种流动类型的递推公式组合起来，可求解单一河道中流向顺逆不定的各种对流输移问题。

2、节点方程的形成及求解

在河网计算中，一般来讲，对于一条不与边界节点相连接的河道，其首、末断面浓度都是未知数，使得一河道上的N个断面浓度未知数，最多只能获得N-1个方程，方程中的未知数多于方程数，无法直接求解，需要利用河道节点方程来补充方程的个数。对于与任意一个河道节点相连的河道首或末断面，如果该断面的流向为流出节点，则取该断面浓度为节点浓度；如果该断面上流向为流入节点，则根据该断面所在河道的递推方程组获得该断面浓度的算式，代入式(31)或式(32)，获得节点浓度方程:

a_i1C₁+a_i2C₂+…+a_imC_m＝r_i (i＝1,2,…,m) (50)

式中m表示概化河网中的节点数，a_ij为节点j的浓度与节点i浓度的相关系数。

一个网河中M个节点可获得M个方程，可形成方程组：

AC＝R (51)

其中：

为系数矩阵；C＝[C₁ C₂ Λ C_m]^T为节点浓度向量，R＝[r₁r₂ Λ r_m]^T为右端常数项。m个节点浓度未知数，m个节点方程，方程组可解。

由于求解的方程组(51)仅涉及到内节点的浓度，而不涉及各断面的浓度，这就大大缩小了系数矩阵A的阶数，相应地节约了计算机内存的需求，并减少了计算时间。

对系数矩阵A，如果节点i与节点j不相邻，或节点i与节点j相连，但连接但与节点i的流量方向并不是流向节点i，则a_ij＝0。对于一般平原河网，与某一内节点i相邻的内节点数，通常是整个河网内节点总数的一小部分，因此节点i的代数方程，仅涉及部分内节点的未知浓度，由此建立的全部内节点的代数方程组，其系数矩阵是一个稀疏矩阵，即矩阵A有多数元素为0。又因为向同一节点汇合的河道只有少数几条，某时刻出现的流动类型也只有几种，这就大大增加了系数矩阵A中的零元素的个数，增大了矩阵的稀疏程度。

为了与水动力模型节点方程组的求解方法保持一致，采用矩阵标识法进行求解。

3、水质模型求解步骤

1)在河网水流计算的基础上，根据河道的流态，建立每条河道上各断面浓度的递推方程组；

2)建立节点浓度方程组；

3)根据节点浓度方程组，求得网河中每个节点的浓度值；

4)将节点浓度值回带给与节点相连的河道首、末断面未知量；

5)利用河道上的递推方程组，求解河道上各断面的浓度值。

在率定南方河网区的流场后，对水质模拟范围的河流进行水质模拟。对模拟范围内的河流进行重新编码，编码与水动力模型的编码具有一一对应关系，故断面流量、水位等水文条件，可以直接采用水动力模型的计算结果。

数据输出文件由水流数学模型及水质模型两部分的输出结果组成，其中潮位和流量结果由水流数学模型输出，污染物浓度结果由水质数学模型输出。

南方感潮河网区水流水质耦合模拟***(SWRNQ)的运行由数据输入和模型库两大模块决定，并以模型库中水流数学模型、水质数学模型两部分的有机结合为基础，最终输出水量或水质结果。水流数学模型对水流模型所需数据进行计算处理，为水质数学模型提供必要的基础数据。

整个***的运行机制过程大致如下：水流数学模型根据输入模块提供的水流模型所需数据文件经过运算，一方面可以输出任意断面的时段水位流量数据，另一方面还可以为水质数学模型提供计算所需的河道断面水位流量过程数据；然后水质数学模型根据输入模块提供的水质模型所需文件，并结合水流数学模型提供的河道断面水位流量数据，在水流水质模型耦合关系的作用下，经运行并输出任意所需断面的污染物浓度变化过程。

需要说明的是：水流数学模型可以独立运行，但水质数学模型则需建立在水流数学模型之上，即由水流数学模型为其提供计算所需的时段河道断面水位流量数据。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (3)

1.南方感潮河网区水流水质耦合模拟***，其特征在于，所述***包括：

数据输入模块：用于输入数据，输入数据由水流模拟、水质模拟和水流水质模型耦合关系三部分所需数据文件组成，其中水流模拟所需文件包括河道流场边界文件，流场边界类型及流场河道节点文件，河道地形数据文件，取水位置和取水量文件；水质模拟所需文件包括浓度场河道首末节点文件，浓度场污染源文件，浓度场外河道地形走向判断因子文件；水流水质耦合关系所需文件为浓度场与流场的河道、节点对应关系文件；

模型库：包括水流数学模型和水质数学模型两大模块，且两者相互连接，有机结合；所述水流数学模型导入和存储河道地形数据、模型边界水文数据、取用水数据和模型河道糙率参数初始值，然后根据输入数据自动进行模型参数率定，最后使用率定后的模型进行工程取水论证，以期达到水资源的可持续开发和利用，并得到水质数学模型所需的基础数据；所述水质数学模型在水流水质模型耦合关系下通过水流数学模型提供的基础数据，结合输入的污染源位置及排污量数据、模型边界水质浓度数据和污染物降解系数等进行河道水质模拟演算，制定合理的排污方案；

数据输出模块：用于输出数据，数据输出文件由水流数学模型及水质模型两部分的输出结果组成，其中潮位和流量结果由水流数学模型输出，污染物浓度结果由水质数学模型输出。

2.根据权利要求1所述的南方感潮河网区水流水质耦合模拟***，其特征在于，所述水流数学模型以一维非恒定流圣维南方程组为基础，采用四点加权Preissmann隐式差分格式对方程进行离散，以提高模型稳定性及模拟精度；采用矩阵标识法对节点方程进行求解，以提高计算效率，实现节点任意编码，增强可扩充性、可移植性和通用性；根据南方感潮河网区河道水闸众多的现状可自行添加水闸，以更真实地反映河网区的水动力过程；对于河道糙率率定提供了一套适用于南方感潮河网区的调整规则。

3.根据权利要求1所述的南方感潮河网区水流水质耦合模拟***，其特征在于，所述水质数学模型采用适用于各种流动情况下求解河网对流输移问题隐式差分方程的一套计算公式及相应的河网计算方法，对于网河区某一河道，根据河道两端的水流方向的不同组合设置了4种流态，即：顺流向的流动、逆流向的流动、河道两端向中间的流动、河道中间任意河段向两端的流动；在河网计算中，采用了一种结合河道节点方程来补充水质方程的个数的求解方法。