一维水量水质模型

1、第七章 一维非恒定河流和河网水量水质模型对于中小型河流， 通常其宽度及水深相对于长度数量较小， 扩散质（污染物 质、热量）很容易在垂向及横向上达到均匀混合， 即扩散质浓度在断面上基本达 到均匀状态。这种情况下， 我们只需要知道扩散质在断面内的平均分配状况， 就 可以把握整个河道的扩散质空间分布特征， 这是我们可以采用一维圣维南方程描 述河流水动力特征或水量特征（水位、流量、槽蓄量等） ；用一维纵向分散方程 描述扩散质在时间及河流纵向上的变化状况。 特别地， 对于稳态水流， 可以采用 常规水动力学方法推算水位、 断面平均流速的沿程变化； 采用分段解析解法计算 扩散质浓度沿纵向的变化特征。 但是，

2、在非稳态情况下 （水流随时间变化或扩散 质源强随时间变化） 解析解法将无能为力 （水流非恒定） 或十分繁琐 （水流稳态、 源强非恒定），这时通常采用数值解法求解河道水量、水质的时间、空间分布。 在模拟方法上， 无论是单一河道还是由众多单一河道构成的河网， 若采用空间一 维手段求解， 描述水流、水质空间分布规律的控制方程是相同的， 只不过在具体 求解方法上有所差异而已。7.1单一河道的控制方程7.1.1水量控制方程采用一维圣维南方程组描述水流的运动，基本控制方程为：q xzb wtq(1)qtq2u gaxz2xuagqnu4/302xr(2)式中t为时间坐标，x为空间坐标，q为断面流量，z为断

3、面平均水位，u为断面平均流速，n为河段的糙率，a为过流断面面积，bw为水面宽度（包括主流宽度 及仅起调蓄作用的附加宽度） ，r为水力半径，q为旁侧入流流量（单位河长上 旁侧入流场）。此方程组属于二元一阶双曲型拟线性方程组，对于非恒定问题， 现阶段尚无法直接求出其解析解， 通常用有限差分法或其它数学离散方法求其数 值解。在水流稳态、 棱柱形河道条件下， 上述控制方程组退化为水力学的谢才公 式，可采用相应的方法求解水流特征。7.1.2扩散质输运控制方程 描述河道扩散物质运动及浓度变化规律的控制方程为： 带源的一维对流分散（弥散）方程，形式如下：（ac）（qc）aexxxcrx hakacsrs（3

4、）t x式中，c为污染物质的断面平均浓度，q为流量， ex为纵向分散系数，s为单位 时间内、单位河长上的污染物质排放量，k为污染物降解系数，sr为河床底泥释 放污染物的速率。1此方程属于一元二阶偏微分方程， 对于非恒定水流问题， 微分方程位变系数 的偏微分方程， 现阶段尚无法直接求出其解析解， 通常用有限差分法或其它数学 离散方法求其数值解。 在水流稳态、 污染源源强恒定条件下， 可按水动力特征将 河道分为若干子段，在每个分段上，上述控制方程简化为常系数的常微分方程， 可采用解析方法秋初起理论解。7.2单一河道一维水量水质模型7.2.1单一河道一维水量模型(1)控制方程的离散采用四点隐式差分格

5、式离散方程组。如图1所示，河道被(n+ 1)个断面分为n个子河段，在第i个子河段m(i,i+ 1)上，对任一变量 取：(m ) (ij(m)ij1) / 2(4)jij 1!ij 1(1 )i 1 i(5)xx(m)tj 1 j 1 j ji i 1 i i 1(6)2t第i个子河段i i+1 n n+1图1计算断面示意图式中，上角标表示时间坐标，下脚标表示空间坐标。 为空间差商的权重系 数(01),=0时，此格式为显式格式， 而当0时，此格式具有隐式差分的特征。为使差分方程保持无条件稳定， 必须0.5。采用下式进行阻力项的线性 化：r r ir i 1nqug4/3g4/3qi2n u2j

6、1n u2j 14/3qi10.5(7)j1j 1将式(4)-(6)代入连续方程得第i个子河段的差分方程：ciziciz i 1qiqi 1di(8)xi(1 )i2t下角标i+1/2表示断面i与断面i+1河段的均值。按照同样的方法，可得动量方 程的差分方程：式中，cibwi 1/2，di(qiqi 1) ci(zijjjzi1)jqixieiqigiqi 1fizifizi 1hi(9)2式中，eifigihi对任一河段xit2ui 1/22gaxig xi2n24/32jbu i 1/ 22t2ui1/2jg xi2ui 1/2n24/ 3i1jx itqi 1/2j21(qijqi1)j

7、12 j jga bu i 1/2(zi 1z1)j2jji(in)，可得到方程组：cizi 1giqiqiqi 1dicizi(10)1fizifizi 1hi对每一河段可列出两个线性代数方程， 再加上上eiqi下游边界条件， 构成完备的封闭 方程组，采用追赶法可求得各个断面的水位流量。2)边界条件 根据上有下游边界条件类型的不同可以写成如下两种追赶形式：(上游水位边界条件 z1z1*(t) ；下游水位(或流量)边界条件zn1zn 1*(t)(或qn 1qn 1(t)，追赶形式为：z1p1r1q1q1l2m2q2z2p2r2q2li 1m i 1qi 1qizi 1pi 1ri 1qi 1q

8、nln 1m n 1qi 1zn 1pn 1rn 1qn 1zn 1zn* 1(t) (or,qn 1qn* 1(t)式中， p,r,l,m 为已知系数，依据上述方程组，可逐步由下边界水位或者流量，推算得到上游各个断面水位流量值*(11)上游流量边界条件 q1q1(t ) ；下游水位边界条件zn 1zn 1(t)，追赶形式为：*3q1p1r1z1z1l2m2z2q2p2li 1r2z2ziqi 1m i 1zi 1pi 1ri 1zi 1(12)znqn 1zn 1ln 1m n 1zi 1pn 1rn 1zn 1zn 1(t)*式中， p,r,l,m 为已知系数，依据上述方程组，可逐步由下边

9、界水位，推算得 到上游各个断面水位流量值。7.2.2单一河道一维水质模型（1）控制方程的离散与求解 对方程（3）进行离散，空间差分采用隐式迎风差分格式。顺流时（从断面i流 向i+1）有：(ac) t(ac)i(ac)itn(auc)(qc)i(qc)i 1xi 1x(ae1( aexi 1ki 1ai 1/2cici 1)ixici(aex)i 1cici 1xxi 1kacsi 1得到统一形式的差分方程：aicibicicici 1di(13)式中， ai,bi,ci,di为系数，分别表示为：biqi1/2,0(aex)i 1/ 2 ti 1/2xi 1aiviqi 1/ 2,0 t qi

10、1/ 2,0 t(aex)i 1/2xi(aex)i 1/ 2txi 1kivitci qi 1/2,0(aex)i 1/2 txi4diviciqi,0 sixit qi,0si 1方程(13)两边同nnxi 1t时除以 ai得到：ci在顺流情况下eici 1fici 1gi，各河段差分方程可写成：(14)c2e2c1f2c3g2c3e3c2f3c4g3cieici 1fici 1gi(15)cn 1en 1cn 2fn 1cngn 1cnencn 1fncn 1gn对首断面给定第一类边界条件， 对末断面给定第二类边界条件，可得到如下封闭的方程组：c1c1(t)*c2c3cie2c1f2c3

11、g2e3c2f3c4g3(16)eici 1fici 1gicn 1en 1cn 2fn 1cngn 1cnencn 1en 1cnfncn 1gn 1gncn 1对方程组（16）采用追赶法可容易求得 n,n 1,n 2, ,3,2 等断面的扩散质的浓度。（2）参数确定纵向分散系数e的确定e与水流流速、水面宽度成正比，与水深成反比，常采用下面xx的经验公式：exc0q式中，c0c/ g是无尺度谢才系数，c为谢才系数，=b/ h为宽深比，q为单 宽流量，0.01为经验常数。降解系数k的确定5可采用监测资料对降解系数进行率定，或根据经验得到7.2.2应用实例三峡大坝位于宜昌县三斗坪中堡岛,葛洲坝位

12、于南津关下游的宜昌市境内.两坝间水域处于鄂西山区向平原的过渡地带,周围地形地貌呈西北高东南低之势 三峡大坝至葛洲坝河段长38km,两坝间江面宽210m至1500m,大部分处于西陵峡 谷中。该实例建立了两坝间水量水质模型,分别运用大坝一期围堰及二期围堰施 工期间的同步水文水质实测资料对模型进行了率定和验证,取得了较好效果.此 模型可用以预测大坝施工期间及投入运行后两坝间水流及水质特性的变化.（）水量模型率定利用1996年三斗坪、白庙子及黄陵庙等断面水文观测资料率定糙率，得到 各子河段的糙率， 率定结果显示糙率分布取值范围为0.036至0.050。1996年白 庙子及黄陵庙等断面实测水位过程线及计

13、算水位过程线图略。 （）水量模型验证采用太平溪断面1998年实测流量作模型验证的上边界条件，采用葛洲坝坝 前断面1998年实测水位作模型验证的下边界条件， 对1998年两坝间水流进行模 拟，通过对各水文观测断面的水文要素的观测值及计算值进行比较可见， 吻合程 度较好，因此此水量模型可用来模拟两坝间的一维水流情况，1998年白庙子及 黄陵庙等断面实测水位过程线及计算水位过程线见图。图 2 模型验证各断面实测及计算水位过程线（） codmn模型验证 高锰酸盐指数的降解系数由实测资料，取经验值，不采用模型率定。采用1998年712月的codmn浓度监测资料对codmn模型参数可靠性进行验证，模 拟时

14、codmn污染源分点源和面源两种情况。 部份断面codmn水质因子实测及计 算浓度值见图3。6三斗坪断面cod(mg/l)3.00实测值计算值cod(mg/l)白庙子断面2.001.000 30 60 90 120 150 180 t(day)cod(mg/l)3.00南津关断面实测值 计算值2.001.000 30 60 90 120 150 180t(day)图 3 模型验证codmn浓度实测及计算值7.3平原河网水量水质模型河网地区是中国社会、 经济、文化的发达地区， 在国家经济发展中占有举足轻重的地位。 随着地区经济的进一步发展、居民物质精神生活水准的进一步提高，水资源问题日益突出。

15、如中国著名的长江三角洲、珠江三角洲地区， 已成为水质型缺水区，水资源、水环境问题已 成为制约经济社会发展、 事关区域可持续发展的重要因素。 生产实践的迫切需要使得河网地 区的水环境保护研究呈不断深入、 系统之势。 作为基础， 河网水流水质模拟方法是进行区域 环境规划、 环境管理等的必备工具， 在水环境问题研究中占有重要地位。因此，该领域的研 究一直是环境科学研究人员、 环境管理决策部门十分关注的重要问题。 在物理模型、 数学模 型这两大手段中，因为数学模型具有经济、快捷、 实用等优点，加之水网地区河道密布等客 观条件的限制， 现阶段只能采用数值方法模拟水网地区的水流运动及污染物输运规律。 但由

16、 于理论、技术及各种客观条件的制约， 目前水流特别是水质数值模拟精度不是十分令人满意。7.3.1平原河网水量模型按河网水流的控制方程及对河网的概化处理方式不同， 河网地区水流数值模拟方法可分为两大类： 第一类为常用的一维圣维南方程组数值解法；第二类为所谓的 “组合单元解法” 。，其中， 一维圣维南方程组数值解法又可分为直接解法和间接解法两种。 在直接解法中， 较有 代表性的有文献g.noseda, mathematical model of unsteady flow in open channels networks, proceedings of the international sym

17、posiumon unsteady flow in open channels,1976，j.j.r. williams, t.r.e. chidley, nonliear analysis of unsteady flow in open channelnetworks, proceedings of the international symposium on unsteady flow in open channels,1976提出的方法， 该类方 法将计算断面交替取为水位和流量断面， 对河网在所有计算断面上统一对一维圣维南方程组 差分离散并求解。 但该方法未知数数量较多， 在河网规模较

18、大的情况下， 因为河道的交叉衔 接，形成的矩阵是一个不规则、不对称的大型稀疏矩阵。为减少存贮， j.j.dronkers 于 1976 年提出间接解法5的思想，以后又有许多学者对其作了进一步完善。为提高计算效率，需 缩小矩阵规模， 中山大学数力系 1976年提出了河网非恒定流隐式方程组稀疏矩阵解法6李岳生 等，河网不恒定流隐式方程组稀疏矩阵解法，中山大学学报（自然科学版） ，1979年3月，该方法从河网矩阵本身的特点出发，能够 有效地节省存储并提高计算速度，但矩阵中需包含所有断面的未知数，方程规模仍然较大， 其实际使用也受到限制。 而间接解法是将断面未知数往交汊点7集中， 待求出交汊点未知数后

19、， 再求解各单一河道未知变量，计算效率较高。间接解法的思想首先由荷兰水力学专家 dronkers 于 1976 年提出，以后又有许多学者对此方法进行进一步的完善，相继提出了河网 非恒定流的二级解法徐正凡，明渠非恒定流m.武汉水利电力学院，1983年.、三级解法8张二骏，河网非恒定流的三级联合解法，华东水利 学院学报，1982年第1卷：112和四级解法9吴寿红，河网非恒定流四级解法，水利学报，1985年，第8期：4250。此类间接解法中， 以三级解法最为常用 。其基本求解思路可概括为 “单一河道交汊点单一河道” ，即（ 1） 先将单一河道划分为若干子河段，在计算断面上对一维 saint-vena

20、nt 方程组进行有限差分运 算，得到各单一河道差分方程组，进行消元计算， 得到单一河道首、末断面间流量与水位的 相互关系；（ 2）根据河道交汊点水量守恒方程， 得到并求解交汊点水位方程组， 得到所有交 汊点水位；（ 3）根据交汊点水位返回各单一河道，求得各计算断面水位、流量值。“组合单元解法”10由法国水力学专家 jean a. cunge 于 1975 年首次提出， 国内也有研 究者采用此方法进行了水网地区的水力模拟1112。此方法的基本思想是：将河网地区水力 特性相似、 水位变幅不大的水体概化成单元。 取单元中心的水位为代表水位， 采用谢才公式 模拟单元间流量交换，根据水量守恒建立每一单元

21、的微分形式的水量守恒方程， 离散并得到 以单元水位为自变量的代数方程，辅以边界条件，可求得各单元水位、单元间流量。在上述两类方法中， 组合单元解法对河道进行了简单概化， 以单元为计算单位， 计算相 对简单，但模拟精度相对较低，仅实用于大尺度水域的水力模拟。而 saint-venat 方程组数 值解法可以精确计算每一条河道的水流状况， 所以成为目前河网水力模拟的主流方法， 其中 又以三级解法最为常用13141516。 当采用三级解法模拟长系列大范围河网水流特性时， 需 对河道进行概化处理，即将等级较小的河道概化为一条“概化”河道，要求该概化河道与被 概化河道的过流能力、调蓄量相当。 同时，为考虑

22、降雨对河网水力特性的影响，还必须对河 道包围的陆域面积进行产汇流计算， 将产流量以包围陆域的河道长度为权重分配到周围河道17。由于三级解法以单一河道为模拟对象，同时考虑了降雨等的影响，计算精度较高，可 以满足水利、航运、环境保护等的需要。（1）控制方程的离散与求解采用四点隐式差分格式离散方程，得子河段差分方程：cizicizi 1qiqi 1di（17）eiqigiqi 1fizifizi 1hi式中ci,di,ei,fi,gi,hi由时段初值及河道特征求得。（18）（2）节点连接条件水流运动在河网各节点上应满足质量守恒及能量守恒，即满足以下两个连接条件：质量守恒条件进出某一节点的流量与该节点

23、内水量蓄量的增减相平衡，定量表示为：mqi（ki1j 1kj） / t（19）式中 k 为节点编号， m 为流入 （流出）节点 k 的河道数，k为节点蓄量， q 为流量（流入为 正；流出为负 )。若节点汇合区容积与子河段容积相比可忽略不计，则此节点称为无调节节m点，方程 (19)可简化为：i1qi 0能量守恒条件不计节点处能量损失有：zizji =1,2,m, j 1,2, m(20)8(3)边界条件有三种类型的边界条件：水位边界条件即在边界河道上给定水位随时间的变化过程：z z(t)；流量边界条件即在边界河道上给定流量随时间的变化过程： q q (t ) ；水位流量关系 当边界河道上有水工建

24、筑物(如水闸、堰、堤坝等)时，通常给定水位流量关系：q q ( z )。4) 方程的求解利用消元法，方程 (3) 、(4)经递推运算，写成如下形式：qii iziiz n 1(21)(22)qi 1 i 1 i 1zi 1 i 1z1式中，i,iiiii为追赶系数，可逐步递推求得。利用公式, , , ,(7)依次由末断面向首断面递推，可将各断面流量表为该断面水位及末断面水位的函数。利用公式 (22) 依次由 首断面向末断面递推， 可将各断面流量表为该断面水位及首断面水位的函数。河道首末断面分别与节点相连，分别对应如下追赶方程：q1特殊的， 单一1 1z1 1z n 1(23)n 1z1(24)

25、qn 1n 1 n 1z n 1将节点各支流相应的(23)或(24)代入公式(19)，并将节点各相邻断面水位统一表示 成节点水位，得节点方程：f (z ,zi i邻1,z邻2,z邻ki) 0i=1,2, ,m式中，m为节点总数，ki为第i个节点相连河道数。当某单一河道为边界河道时， 分三种类型分别有： 对于第一类边界条件， 式(24)中 z1已知；对于第二类边界条件，式(23)中 q1已知；对于第三类边界条件，已知线性关系 q1=q(z1)，代入(23)消去 q1。无论何种边界条件，最终可增加一9水位方程，封闭节点方程组fi（zi,z邻1,z邻2,z邻ki）位，再返回单一河道方程，最终求得各断

26、面水位及流量。0，可解得各节点水7.3.2平原河网水质模型（1）单一河道控制方程的离散采用如下网格对方程（3）进行离散，空间差分采用隐式迎风差分格式，由于 河网地区水流流向可能谁随时间变化， 根据流向的不同， 分顺流及逆流两种状况 给出差分形式：ci 1cici 1xi 1xi图4水质模拟网格布置示意图顺流时（从断面i流向i+1）有：(ac)(ac)i(ac)intt(auc)(qc)i(qc)i 1xxi 1c1ci 1cixic c)i 1i i 1xi 1(ae)x)ii i 1(aexxx( aexxi 1kacskis1ai 1/2cii 1逆流时（从断面i+1流向i）有：(ac)(

27、ac)i(ac)intt(auc)(qc)i1(qc)ixxicci 1ci(aexxi 1x)1xx( aex(ae)i 1cici 1x)iixikacs kiai1/2cisi为考虑流向顺逆变化的影响，引入流向调节因子rc及rd，对于断面表示为：i10qw(qiqi 1) /2qe(qiqi 1) /2(qwqw) /2qwrce(qeqe)/2q (qerdw(qwqw) / 2qwrdeqe)/2qerd0,当qw,qe0时)得到任意流向下，统一形式的差分方程：ici1iciici 1zi式中，i,i,i,zi为系数，对于一般断面(i=2,3,n-1)分别表示为：rcwdwwrdwd

28、pwfcw t /vrcwdwwrd wdpwrcedpprdedepfcpfdpt /vrcwki 1rdekit 1.0ir dce ppfde t /vrd edepzi11circwsi 1rdesi ti 1/2i 1/2对于首断面(i =1)，逆流时有：0rd edeptfdpx ardekrdexiai 1/2ittrdedepfdexiai 1/2tzc11rdes1ai 1/2对于末断面(i =n)，顺流时有：f t fcwrcwdan 1/2wwtxrcwdwwficprcwki 11ai1/2sn 1tzncnnai 1/2式中，ai 1/2(ai1ai) /2 ,ai

29、1/2(ai1ai) /2(25)11v1ai 1/2xi 1,v2ai 1/2xivrcwv1rdev2dww(aex)i1/(aex)i/xi 1,xi 1,)/2,dpp(aex)(aex)ii/xidpwdep/ xifcw(qi 1(qiqi 1fc p(qiqi)/2fd pqi 1) / 2 ,fde(qi 1qi 1) /2将差分方程(25)用向量形式表示为：a cc cbc(26)式中 ac为一(n 1) n阶的三对角的系数矩阵， cc为此河道断面平均浓度的n维列向量，bc为已知的(n-1)维列向量。(26)式中含(n-1)个方程，n个未知数，方 程组不闭合，需引入节点方程及

30、边界条件。(2)河网节点方程对于污染物充分混合的节点，从时刻n t到(n 1) t ，根据节点有无调蓄作 用，分别给出相应的节点方程。若节点有m条单一河道，其中流入节点的河道 m1条，流出节点的河道 m2条，m1及 m2随流场变化而变化(m1+m2=m)。若节点本身具有调蓄作用，有：c out , icn，i=m11,m12, ,m(27)m1qin ,ici1m in ,iqout , icout ,ii m1 1(28)tncnncne sn式中，cn为节点的浓度， cout ,i为流出节点的第i条河道与该节点相邻断面的污 染物浓度， ci n,i为流入节点的第i条河道与该节点相邻断面的污染物浓度，q分别为与之相应的流量，n为节点的蓄量，sn为节点的污染源加入项，方程(28)右端括号中末项为时段初节点污染物降解至时段末的残留量。(27)式称为充分混 合假定，(28)式反映节点内污染物的质量守恒。当节点调蓄可忽略时，称为无调蓄节点，有：12n0(3)模型的动态耦合数值计算在流场已知时，由(26)-(28)式加上边界条件及初始条件，构成封闭的线性代 数方程组，即可求解。但对于大型复杂河网来说，其规模巨大，较为经济的做法 是：将方程(26)作适当的运算，得单一河道首末断面浓度关系方程，并与充分混 合假定(27)及边界条件，一并代入方程组(28)，形成封闭