一维水量水质模型

第七章一维非恒定河流和河网水量水质模型对于中小型河流，通常其宽度及水深相对于长度数量较小，扩散质（污染物质、热量）很容易在垂向及横向上达到均匀混合，即扩散质浓度在断面上基本达到均匀状态。这种情况下，我们只需要知道扩散质在断面内的平均分配状况，就可以把握整个河道的扩散质空间分布特征，这是我们可以采用一维圣维南方程描述河流水动力特征或水量特征（水位、流量、槽蓄量等）；用一维纵向分散方程描述扩散质在时间及河流纵向上的变化状况。特别地，对于稳态水流，可以采用常规水动力学方法推算水位、断面平均流速的沿程变化；采用分段解析解法计算扩散质浓度沿纵向的变化特征。但是，在非稳态情况下（水流随时间变化或扩散质源强随时间变化）解析解法将无能为力（水流非恒定）或十分繁琐（水流稳态、源强非恒定），这时通常采用数值解法求解河道水量、水质的时间、空间分布。在模拟方法上，无论是单一河道还是由众多单一河道构成的河网，若采用空间一维手段求解，描述水流、水质空间分布规律的控制方程是相同的，只不过在具体求解方法上有所差异而已。单一河道的控制方程水量控制方程采用一维圣维南方程组描述水流的运动，基本控制方程为(1)(2)晋+2u磐+gA聖-u2聖+gnRUQ=0at ax axax R4/3(2)式中t为时间坐标，x为空间坐标，Q为断面流量，Z为断面平均水位，u为断面平均流速，n为河段的糙率，A为过流断面面积，BW为水面宽度（包括主流宽度及仅起调蓄作用的附加宽度），R为水力半径，q为旁侧入流流量（单位河长上旁侧入流场）。此方程组属于二元一阶双曲型拟线性方程组，对于非恒定问题，现阶段尚无法直接求出其解析解，通常用有限差分法或其它数学离散方法求其数值解。在水流稳态、棱柱形河道条件下，上述控制方程组退化为水力学的谢才公式，可采用相应的方法求解水流特征。扩散质输运控制方程描述河道扩散物质运动及浓度变化规律的控制方程为：带源的一维对流分散（弥散）方程，形式如下：a(AC)aa(AC)a(QC)aat+axA-KAC+—S+Shr(3)式中，C为污染物质的断面平均浓度，Q为流量，E为纵向分散系数，S为单位x时间内、单位河长上的污染物质排放量，K为污染物降解系数，Sr为河床底泥释放污染物的速率。

此方程属于一元二阶偏微分方程，对于非恒定水流问题，微分方程位变系数的偏微分方程，现阶段尚无法直接求出其解析解，通常用有限差分法或其它数学离散方法求其数值解。在水流稳态、污染源源强恒定条件下，可按水动力特征将河道分为若干子段，在每个分段上，上述控制方程简化为常系数的常微分方程，可采用解析方法秋初起理论解。单一河道一维水量水质模型7.2.1单一河道一维水量模型控制方程的离散采用四点隐式差分格式离散方程组。如图1所示，河道被(n+1)个断面分为n个子河段，在第i个子河段M(i,i+1)上，对任一变量g取：g(M)二(gj+gj)/2TOC\o"1-5"\h\z\o"CurrentDocument"i+1ag(m)e(gj+!-gj+J+(1-e)(gj-gJnn+1Q)i+1i i+1 inn+1Q)Ax淹(M) gj+1 +gj+1-gj - gji i+1i i+1’ 2At第i个子河段I’ii+1图1计算断面示意图式中，上角标表示时间坐标，下脚标表示空间坐标。e为空间差商的权重系数(0<e<i),e=0时，此格式为显式格式，而当e丰0时，此格式具有隐式差分的特征。为使差分方程保持无条件稳定，必须e>05。采用下式进行阻力项的线性化：n2Q|u'R4/3n2Q|u'R4/3IR4/3Qj+1+iIR4/3丿.1、jQj+1-0.5i+1丿将式(4)-(6)代入连续方程得第i个子河段的差分方程：(8)CZ+CZ-Q+Q二D(8)ii i i+1 i i+1 iAx (1-e) Ax式中，q=BY"站，£=存(灯-蹄)+9(Z-j+z+1)+q,it下角标i+1/2表示断面i与断面i+1河段的均值。按照同样的方法，可得动量方程的差分方程：EQ+GQ—FZ+FZ二Hii ii+1 iiii+1 i

式中，厂 Ax cE二—一2uji 20式中，厂 Ax cE二—一2uji 20At i+1/2gAx20(a-Bu2)ji+1/2Ax—+2uj+20At i+1/2gAxi20n2uji+1h=0ArQi+i/2+21-0a—j-N-Bu2)+i/2(乞i-勺)对任一河段i(i=1,2, n)，可得到方程组：CZ+CZ-Q+Q=Dii i i+1 i i+1 i (10)EQ+GQ-FZ+FZ=Hii i i+1 i ii i+1 i对每一河段可列出两个线性代数方程，再加上上下游边界条件，构成完备的封闭方程组，采用追赶法可求得各个断面的水位流量。边界条件根据上有下游边界条件类型的不同可以写成如下两种追赶形式：•上游水位边界条件Z二Z*(t)；下游水位(或流量)边界条件Z=Z*(t)(或1 1 n+1 n+1Q=Q*(t))，追赶形式为：n+1 n+1'Z=P+RQiii1Q=L+MQ1222Z=P+RQ.2222Q=L+MQVi i+1 i+1i+1 (11)Z=P+RQ.i+1 i+1 i+1i+1Q=L+MQn n+1 n+1i+1Z=P+RQn+1n+1 n+1n+1Z=Z*(t)(or,Q=Q*(t))n+1 n+1 n+1 n+1式中，P,R,L,M为已知系数，依据上述方程组，可逐步由下边界水位或者流量,推算得到上游各个断面水位流量值。°上游流量边界条件Q=Q*(t)；下游水位边界条件Z=Z*(t)，追赶形式为:1 1 n+1 n+1

Q=P+RZiiiiZ=L+MZ1222Q=P+RZ.2222Z=L+MZvii+1 i+1i+1 仃2'Q=P+RZi+1i+1i+1i+1Z=L+MZnn+1 n+1i+1Q=P+RZn+1n+1n+1n+1Z=Z*(t)n+1n+1式中，P,R,L,M为已知系数，依据上述方程组，可逐步由下边界水位，推算得到上游各个断面水位流量值。单一河道一维水质模型控制方程的离散与求解对方程(3)进行离散，空间差分采用隐式迎风差分格式。顺流时(从断面i流向i+1)有：Q(AC) (AC)-(AC)nTOC\o"1-5"\h\z= i i—dt At°(AUC) (QC)-(QC)= i i1dx Axi-15dC1 C-C C-C\o"CurrentDocument"(AE—)- [(AE)乜4—(AE) *5x x5xAx xiAx xi-1Ax\o"CurrentDocument"i-1 i i-1-KAC+S—-KAC+Si-1i-1/2i i-1得到统一形式的差分方程：aC=bCaC=bCiiii-1+cC+dii+1 i(13)式中，—V+[Q ,0]式中，—V+[Q ,0]At+[—Q ,0]At+i i+1/2 i-1/2(AEx)i++/Axi(AEx)•At+Axi-1^-+^2•At+KVAta,b,c,d为系数，分别表示为：iiii—{[QJ+哙十}心i-1c—{—[Q ,0]+ z2}・Ati i+1/2 Axd=Vn.Cn+[Q,0]•S-AxAt+[—Q,0]S-AxAtTOC\o"1-5"\h\ziii i ii i i—1 i—1方程(13)两边同时除以a得到：i(14)C=EC+FC+G(14)i ii—1 ii+1 i在顺流情况下，各河段差分方程可写成：'C=EC+FC+G21 23 2C=EC+FC+G(15)32 34 3(15)C=EC+FC+Gii—1ii+1 iC=EC+FC+Gn—1 n—1n—2 n—1n n—1C=EC+FC+Gnnn—1 nn+1 n对首断面给定第一类边界条件，对末断面给定第二类边界条件，可得到如下封闭的方程组：C=C*(t)11C=EC+FC+GTOC\o"1-5"\h\z21 23 2C=EC+FC+G3.3 2 3 4 3C=EC+FC+G (16)iii—1ii+1 iC=EC+FC+Gn—1 n—1n—2 n—1n n—1C=EC+FC+Gnnn—1 nn+1 nC=EC+Gn+1 n+1n n+1对方程组(16)采用追赶法可容易求得n,n—1,n—2,…,3,2等断面的扩散质的浓度。参数确定纵向分散系数E的确定XE与水流流速、水面宽度成正比，与水深成反比，常采用下面的经验公式：XE=aC02q

x0式中，C=cI品是无尺度谢才系数，c为谢才系数，0=B/h为宽深比，q为单宽流量，a=0.011为经验常数。降解系数K的确定

可采用监测资料对降解系数进行率定，或根据经验得到。应用实例［］三峡大坝位于宜昌县三斗坪中堡岛,葛洲坝位于南津关下游的宜昌市境内.两坝间水域处于鄂西山区向平原的过渡地带,周围地形地貌呈西北高东南低之势.三峡大坝至葛洲坝河段长38Km,两坝间江面宽210m至1500m,大部分处于西陵峡谷中。该实例建立了两坝间水量水质模型,分别运用大坝一期围堰及二期围堰施工期间的同步水文水质实测资料对模型进行了率定和验证,取得了较好效果.此模型可用以预测大坝施工期间及投入运行后两坝间水流及水质特性的变化.水量模型率定利用1996年三斗坪、白庙子及黄陵庙等断面水文观测资料率定糙率，得到各子河段的糙率，率定结果显示糙率分布取值范围为0.036至0.050。1996年白庙子及黄陵庙等断面实测水位过程线及计算水位过程线图略。水量模型验证H(m)H(m) 计算值64.00黄陵庙断面72.00"70.00-68.00一66.00"60120180240300360T(daY)0实测值 计算值74.00-72.00-70.00-68.00-66.00-64.000三斗坪断面实测值采用太平溪断面1998年实测流量作模型验证的上边界条件，采用葛洲坝坝前断面1998年实测水位作模型验证的下边界条件，对1998年两坝间水流进行模拟，通过对各水文观测断面的水文要素的观测值及计算值进行比较可见，吻合程度较好，因此此水量模型可用来模拟两坝间的一维水流情况，1998H(m)H(m) 计算值64.00黄陵庙断面72.00"70.00-68.00一66.00"60120180240300360T(daY)0实测值 计算值74.00-72.00-70.00-68.00-66.00-64.000三斗坪断面实测值60 120 180 240 300 360T(daY)图2模型验证各断面实测及计算水位过程线codm模型验证高锰酸盐指数的降解系数由实测资料，取经验值，不采用模型率定。采用1998年7〜12月的CODMn浓度监测资料对CODMn模型参数可靠性进行验证，模拟时CODMn污染源分点源和面源两种情况。部份断面CODMn水质因子实测及计算浓度值见图3。

三斗坪断面COD(mg/L)黄陵庙断面图3模型验证COD三斗坪断面COD(mg/L)黄陵庙断面图3模型验证COD7.3平原河网水量水质模型河网地区是中国社会、经济、文化的发达地区，在国家经济发展中占有举足轻重的地位。随着地区经济的进一步发展、居民物质精神生活水准的进一步提高，水资源问题日益突出。如中国著名的长江三角洲、珠江三角洲地区，已成为水质型缺水区，水资源、水环境问题已成为制约经济社会发展、事关区域可持续发展的重要因素。生产实践的迫切需要使得河网地区的水环境保护研究呈不断深入、系统之势。作为基础，河网水流水质模拟方法是进行区域环境规划、环境管理等的必备工具，在水环境问题研究中占有重要地位。因此，该领域的研究一直是环境科学研究人员、环境管理决策部门十分关注的重要问题。在物理模型、数学模型这两大手段中，因为数学模型具有经济、快捷、实用等优点，加之水网地区河道密布等客观条件的限制，现阶段只能采用数值方法模拟水网地区的水流运动及污染物输运规律。但由于理论、技术及各种客观条件的制约，目前水流特别是水质数值模拟精度不是十分令人满意。7.3.1平原河网水量模型按河网水流的控制方程及对河网的概化处理方式不同，河网地区水流数值模拟方法可分为两大类：第一类为常用的一维圣维南方程组数值解法；第二类为所谓的“组合单元解法”。，其中，一维圣维南方程组数值解法又可分为直接解法和间接解法两种。在直接解法中，较有代表性的有文献［G.Noseda,MathematicalModelofUnsteadyFlowinOpenChannelsNetworks,ProceedingsoftheInternationalSymposiumonUnsteadyFlowinOpenChannels,1976］，［J.J.R.Williams,T.R.E.Chidley,NonliearAnalysisofUnsteadyFlowinOpenChannelNetworks,ProceedingsoftheInternationalSymposiumonUnsteadyFlowinOpenChannels,1976］提出的方法，该类方法将计算断面交替取为水位和流量断面，对河网在所有计算断面上统一对一维圣维南方程组差分离散并求解。但该方法未知数数量较多，在河网规模较大的情况下，因为河道的交叉衔接，形成的矩阵是一个不规则、不对称的大型稀疏矩阵。为减少存贮，JJ.Dronkers于1976年提出间接解法［5］的思想，以后又有许多学者对其作了进一步完善。为提高计算效率，需缩小矩阵规模，中山大学数力系1976年提出了河网非恒定流隐式方程组稀疏矩阵解法［6李岳生等，河网不恒定流隐式方程组稀疏矩阵解法，中山大学学报（自然科学版），1979年3月］，该方法从河网矩阵本身的特点出发，能够有效地节省存储并提高计算速度，但矩阵中需包含所有断面的未知数，方程规模仍然较大，其实际使用也受到限制。而间接解法是将断面未知数往交汊点集中，待求出交汊点未知数后，再求解各单一河道未知变量，计算效率较高。间接解法的思想首先由荷兰水力学专家Dronkers于1976年提出，以后又有许多学者对此方法进行进一步的完善，相继提出了河网非恒定流的二级解法［ ］、三级解法［8卄1一疋丫丿山口丿」族陆丫厶［徐正凡，明渠非恒定流［M］•武汉水利电力学院，1983年八一族陆丫厶［张二骏，河网非恒定流的三级联合解法，华东水利］和四级解法［9 ］。此类间接解法中，以学院学报，1982年第1卷：1~12 吴寿红，河网非恒定流四级解法，水利学报，1985年，第8期：42~50三级解法最为常用［］［］［］。其基本求解思路可概括为“单一河道—交汊点—单一河道”，即（1）先将单一河道划分为若干子河段，在计算断面上对一维Saint-Venant方程组进行有限差分运算，得到各单一河道差分方程组，进行消元计算，得到单一河道首、末断面间流量与水位的相互关系；（2）根据河道交汊点水量守恒方程，得到并求解交汊点水位方程组，得到所有交汊点水位；（3）根据交汊点水位返回各单一河道，求得各计算断面水位、流量值。“组合单元解法”［10］由法国水力学专家JeanA.Cunge于1975年首次提出，国内也有研究者采用此方法进行了水网地区的水力模拟［11］［12］。此方法的基本思想是：将河网地区水力特性相似、水位变幅不大的水体概化成单元。取单元中心的水位为代表水位，采用谢才公式模拟单元间流量交换，根据水量守恒建立每一单元的微分形式的水量守恒方程，离散并得到以单元水位为自变量的代数方程，辅以边界条件，可求得各单元水位、单元间流量。在上述两类方法中，组合单元解法对河道进行了简单概化，以单元为计算单位，计算相对简单，但模拟精度相对较低，仅实用于大尺度水域的水力模拟。而Saint-Venat方程组数值解法可以精确计算每一条河道的水流状况，所以成为目前河网水力模拟的主流方法，其中又以三级解法最为常用［13］［14］［15］［16］。当采用三级解法模拟长系列大范围河网水流特性时，需对河道进行概化处理，即将等级较小的河道概化为一条“概化”河道，要求该概化河道与被概化河道的过流能力、调蓄量相当。同时，为考虑降雨对河网水力特性的影响，还必须对河道包围的陆域面积进行产汇流计算，将产流量以包围陆域的河道长度为权重分配到周围河道［17］。由于三级解法以单一河道为模拟对象，同时考虑了降雨等的影响，计算精度较高，可以满足水利、航运、环境保护等的需要。（1）控制方程的离散与求解采用四点隐式差分格式离散方程，得子河段差分方程：TOC\o"1-5"\h\zCZ+CZ—Q+Q=D （17）ii ii+1 i i+1 iEQ+GQ—FZ+FZ=H （18）ii ii+1 iiii+1 i式中C,D,E,F,G,H由时段初值及河道特征求得。iiiiii（2）节点连接条件水流运动在河网各节点上应满足质量守恒及能量守恒，即满足以下两个连接条件：•质量守恒条件进出某一节点的流量与该节点内水量蓄量的增减相平衡，定量表示为：£Q—（0j+1—0j）/At （19）iKKi—1式中K为节点编号，m为流入（流出）节点K的河道数，0为节点蓄量，Q为流量（流入为K正；流出为负)。若节点汇合区容积与子河段容积相比可忽略不计，则此节点称为无调节节点，方程(19)可简化为：Q=0ii=1•能量守恒条件不计节点处能量损失有：Z二Z i=i,2,…m, j二1,2，•…m (20)ij边界条件有三种类型的边界条件：水位边界条件即在边界河道上给定水位随时间的变化过程：Z二Z(t)；流量边界条件即在边界河道上给定流量随时间的变化过程：Q=Q(t)；水位流量关系当边界河道上有水工建筑物(如水闸、堰、堤坝等)时，通常给定水位流量关系：Q=Q(Z)。方程的求解利用消元法，方程(3)、(4)经递推运算，写成如下形式：TOC\o"1-5"\h\zQ=a+卩Z+YZ (21)i i ii in+1Q二g+Cz+nz (22)i+1 i+1 i+1i+1 i+11式中，a,卩，丫,g，匚,n为追赶系数，可逐步递推求得。利用公式(7)依次由末断面iiiiii向首断面递推，可将各断面流量表为该断面水位及末断面水位的函数。利用公式(22)依次由首断面向末断面递推，可将各断面流量表为该断面水位及首断面水位的函数。特殊的，单一河道首末断面分别与节点相连，分别对应如下追赶方程：TOC\o"1-5"\h\zQ=a+卩Z+丫Z (23)1 1 11 1n+1Q-C+匚Z+nZ (24)n+1 n+1 n+1n+1 n+11将节点各支流相应的(23)或(24)代入公式(19)，并将节点各相邻断面水位统一表示成节点水位，得节点方程：f(Z，ZZ…Z)-0 i=1，2，…，Mii邻1,邻2, 邻KI式中，M为节点总数，KI为第i个节点相连河道数。当某单一河道为边界河道时，分三种类型分别有：对于第一类边界条件，式(24)中Z已知；对于第二类边界条件，式(23)中Q已知；对于第三类边界条件，11已知线性关系Q=Q(Z)，代入(23)消去Q。无论何种边界条件，最终可增加一111

水位方程，封闭节点方程组f(Z,ZZZ)—0,可解得各节点水邻1,邻2, 邻KI位，再返回单一河道方程，最终求得各断面水位及流量。平原河网水质模型(1)单一河道控制方程的离散采用如下网格对方程(3)进行离散，空间差分采用隐式迎风差分格式，由于河网地区水流流向可能谁随时间变化，根据流向的不同，分顺流及逆流两种状况给出差分形式：Ci+Ci+1i+1n-1 ni-13…i-1TOC\o"1-5"\h\zAx Axi-1 i图4水质模拟网格布置示意图顺流时(从断面i流向i+1)有：d(AC) (AC)-(AC)n\o"CurrentDocument"= i i—dt At°(AUC) (QC)-(QC) = i i~1dx Axi-1d dC1 C-C C-C(AE )- [(AE)乜4-(AE)#Tdx xoxAx xiAx xi-1Axi-1 i i-1—KAC+S=-KAC+Si-1i-1/2i i-1逆流时(从断面i+1流向i)有o(AC) (AC)-(AC)n\o"CurrentDocument"= i i—dt Atd(AUC) (QC) -(QC)= i+1 idx AxiddC1 C-C C-C\o"CurrentDocument"(AE )- [(AE)七4-(AE) 4dxxdx Ax xiAx xi-1Axi i i-1—KAC+S=-KA C+Sii+1/2ii为考虑流向顺逆变化的影响，引入流向调节因子r及r，对于断面i表示为:cdIQ=(Q+Q)/2TOC\o"1-5"\h\zw i i-1r=(Q+Q|)/2QJcw w w wr=(Q-\Q|)/2Qdw w w wr=r=0,cd得到任意流向下，统一形式Q=(Q+Q)/2TOC\o"1-5"\h\ze i i+1r=(Q+|Q|)/2Q

ce e e er=(Q-|Q|)/2Qde e e e(当Q,Q=0时)weaC+卩C+yCii-1 iiii+1=Zi(25)式中，a,卩，丫，Z为系数，对于一般断面(i=2,3,…,n-1)分别表示为:iiiia=-[r D +r D +F]-At/VicwwwdwpwcwB=-[r D +r D +rD+rD +F一F]-At/Vicwwwdwpwcepp deepcpdp+rK+rK]•At+10TOC\o"1-5"\h\zcwi-1 deiy=-[rD+rD-F]-At/Vicepp deep deZ=Z=-Cn+[rSiicwi-1Ai+1/2+rS・ ]Ai+1/2A deii-1/2对于首断面(i=1)，逆流时有:Ia=01B=[rD1deepB=[rD1deepy=-[rd1deepZ=Cn+rS11At-F] +rdpAxAde

ii+1/2

At-F]—deAxAii+1/2Atde1Ai+1/2对于末断面(i=n)，顺流时有：Ata=-[rD+F] ncwwwcwAn-1/2c AtP=[rD+F] +rKAt+rncwww cpAxA cwi-1 cwi-1i-1/2y=0nAtZ=Cn+rS n ncwn-1Ai-1/2式中，A=(A+A)/2i+1/2 i+1 iA=A=(A+A)/2i+1/2 i+1 iTOC\o"1-5"\h\zV=A -Ax1 i-1/2 i-1V=A -Ax2 i+1/2 iV二rV+rVcw1de2TOC\o"1-5"\h\z'D=(AE) /Axww xi-1 i-1D=(AE)/Axpw xi i-1F=(Q+Q|)/2,cw i-1 i-1F=(Q-Q|)/2,Jdp i i-1D=(AE)/Axpp xiiD=(AE) /Axep xi+1 iF=(Q+Q|)/2cp i iF=(Q-Q|)/2de i+1 i+1AC=B (26)cc c式中A为一Ln-1)xn］阶的三对角的系数矩阵，C为此河道断面平均浓度的ncc维列向量，B为已知的(n-1)维列向量。(26)式中含(n-1)个方程，n个未知数，方c程组不闭合，需引入节点方程及边界条件。(2)河网节点方程对于污染物充分混合的节点，从时刻nAt到(n+1)At，根据节点有无调蓄作用，分别给出相应的节点方程。若节点有m条单一河道，其中流入节点的河道m条，流出节点的河道m条,12m及m随流场变化而变化(m+m=m)。1212若节点本身具有调蓄作用，有：C =C， i=m+1,m+2, ,mout,i N 1 1(27)瓦QC-込QCiniini outioutii=1 i=m1+1 (28)GC-QjCje-kat-SAt NNNN N式中，CN为节点的浓度，C为流出节点的第i条河道与该节点相邻断面的污N out,i染物浓度，C为流入节点的第i条河道与该节点相邻断面的污染物浓度，Q分in,i别为与之相应的流量，Q为节点的蓄量，SN为节点的污染源加入项，方程(28)NN右端括号中末项为时段初节点污染物降解至时段末的残留量。(27)式称为充分混合假定，