CN101956382B - 用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法 - Google Patents

Abstract

用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法，是将圣维南基本微分方程组在Δt时段简化为水量平衡方程式和动力方程式后，根据水量平衡方程式和动力方程式分别计算河段槽增量：一是根据水量平衡方程式可计算槽蓄增量，二是根据时段初、末的两个动力方程式也可计算槽蓄增量，由于两种方法计算的河段槽蓄增量应相等，因此本发明就为水力学方法求解圣维南方程组模拟洪水演进提供了一种新的、可靠的方法，在河道、渠道、分蓄洪区和水库合理确定防洪工程规模、洪水预报、科学调度洪水等方面发挥着重要作用。

Description

用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法

技术领域

本发明涉及防洪减灾科学技术领域，特别是根据水量平衡方程式与动力方程式的槽蓄增量关系求解圣维南方程组模拟洪水演进的方法。

背景技术

我国是一个洪涝灾害频繁发生的国家，全国1/2的人口，1/3的耕地和主要大城市都处于江河洪水位以下，受洪水威胁的地区工农业产值占全国2/3。为了达到减轻洪水灾害以及利用水资源为国民经济与社会发展服务的目的，建立模拟洪水演进的数学模型，在河道、渠道、分蓄洪区和水库合理确定防洪工程规模、洪水预报、科学调度洪水等方面发挥着重要作用。

描述河渠非恒定流运动规律的一阶拟线性双曲型偏微分方程，在1871年首先由法国学者圣维南建立。在一般情况下，无法求出其普遍的解析解。目前，只能使用近似的计算方法求解。近似的计算方法大致可分为两类：

1、水力学方法：

直接差分法就是其中的一种方法。直接差分法是用偏差商代替偏微商，将基本微分方程离散化为差分方程，求在自变量域x～t平面差分网格上各节点近似数值解的方法。基本微分方程中的偏微商可用不同形式的偏差商即差分格式代替，从而得到不同的差分方程。但从目前的技术手段而言，研究成果一般用于规划设计阶段，在流域水文预报模型中尚未采用。

2、水文学方法：

水文学方法如著名的“马斯京根法”，是将基本微分方程组分别简化为水量平衡方程式和近似为蓄泄方程式，用差分形式合解水量平衡方程式与蓄泄方程式。缺点是：①概化相对较粗，难于详细分析河道内复杂水流情势以及主要特征；②难于提供河段内水位的时空分布和动态模拟洪水在河段内的演进过程。

发明内容

本发明的用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法，是将圣维南基本微分方程组在Δt时段简化为水量平衡方程式和动力方程式后，根据水量平衡方程式和动力方程式可以分别计算河段槽蓄增量：一是根据水量平衡方程式可计算槽蓄增量，二是根据时段初和末的两个动力方程式也可计算槽蓄增量。由于两种方法计算的河段槽蓄增量应相等，因此，这一发明就为水力学方法求解圣维南方程组模拟洪水演进提供了一种新的、可靠的方法。

本发明的用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法，于2009年11月下旬已在湖北省水利水电勘测设计院进行的杜家台分蓄洪区可行性研究报告修编工作中得到了初步应用，通过模拟洪水演进计算，结果发现没有必要兴建原报告的支洪道项目，若实施研究报告，可节省投资2亿余元。

本发明的用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法，属于水力学方法范畴。本方法既与水文学方法和现行的水力学方法有相似的地方，也有不同的地方。

与水文学方法相似的地方是都涉及到了河段槽蓄内容；不同的是水文学方法如“马斯京根法”是将动力方程式近似变为蓄泄方程式后与水量平衡方程式联解，而本方法是直接将水量平衡方程式与动力方程式联解。

与现行的水力学方法相比，相同之处是都涉及到了需要河道地形资料和依据水文资料率定河段糙率；不同的是，本方法将圣维南基本微分方程组在Δt时段简化为水量平衡方程式和动力方程式后，以河段槽蓄增量的关系联解方程组。

具体方法如下：

现将圣维南基本微分方程组在Δt时段简化为水量平衡方程式(1)和动力方程式(2)以及动力方程式(3)：

${\mathrm{\ΔZ}}_1=\frac{{\stackrel{\‾}{Q_1}}^2}{{\stackrel{\‾}{K_1}}^2}L---\left(2\right)$

${\mathrm{\ΔZ}}_2=\frac{{\stackrel{\‾}{Q_2}}^2}{{\stackrel{\‾}{K_2}}^2}L---\left(3\right)$

以上各式中：L为河段长度，单位m；Δt为时段，单位s；Z为水位，单位m；Q为流量，单位m³/s；W为河段槽蓄水量，单位m³；ΔZ为河段上、下游计算断面的水位差，单位m；

为河段平均流量；ΔW＝(W2-W1)；为河段平均流量模数；各量脚注1、2分别为时段初、末的时刻。

为了说明联解式(1)和式(2)以及式(3)的方法，根据上述方程组绘制成河段槽蓄增量ΔW水量平衡和ΔW动力的图形。

如图2所示是根据水量平衡方程式(1)绘制的，槽蓄增量ΔW水量平衡是在时段Δt内河段总入流量与总出流量的差值，即由入流过程Q_上1～Q_上2与出流过程Q_下1～Q_下2所包围的水体。

如图3所示是根据动力方程式(2)和(3)绘制的，槽蓄增量ΔW动力是在时刻t₁的水面线Z_上1～Z_下1和时刻t₂的水面线Z_上2～Z_下2所包围的水体，即由上游和下游计算断面分别在时刻t₁和t₂的水位转换成过水断面积的差值平均后乘以河段长度得出的体积。

当已知条件为Q_上1、Z_上1、Q_下1、Z_下1、Q_上2以及Q_下2或已知条件为Q_上1、Z_上1、Q_下1、Z_下1、Q_上2以及Z_下2时，可依据槽蓄增量ΔW的关系求解Z_上2和Z_下2或求解Z_上2和Q_下2：

根据方程式(1)，如假定Q_下2或Z_下2值，可以计算出相应的河段槽蓄增量ΔW；

根据方程式(2)、(3)，继续利用上述假定的Q_下2或Z_下2值，也可以计算出河段槽蓄增量ΔW。

通过计算机程序按上述关系调节Q_下2或Z_下2值，使两种方式计算的河段槽蓄增量ΔW基本相等，此时的Z_上2和Z_下2或Z_上2和Q_下2即为方程组的解。

当计算断面个数为n时，则需要增加n-2个断面流量计算关系式求解方程组。

如以上游起始断面的编号为i＝1，依次往下游的断面编号则为i＝2、i＝3…i＝n-2、i＝n-1、i＝n。

当断面编号1＜i＜n时，根据水量平衡方程式(1)，如断面的流量和河段槽蓄增量用二维数组表示，则断面i-1与断面i之间和断面i与断面i+1之间的水量平衡关系式分别为：

Q(i-1，1)+Q(i-1，2)-Q(i，1)-Q(i，2)＝2ΔW(i-1，i)/Δt---------------(4)

Q(i，1)+Q(i，2)-Q(i+1，1)-Q(i+1，2)＝2ΔW(i，i+1)/Δt---------------(5)

联解方程式(4)和(5)：

Q(i，2)＝(ΔW(i，i+1)-ΔW(i-1，i))/Δt+(Q(i-1，2)+Q(i+1，2))/2+(Q(i-1，1)+Q(i+1，1)-2Q(i，1))/2------------------------------(6)

也可近似写成式(7)等形式：

Q(i，2)＝(ΔW(i，i+1)-ΔW(i-1，i))/2/Δt+(Q(i-1，2)+Q(i+1，2))/2------(7)

式中：二维数组流量Q括号内的第一项为断面编号，第二项为时段Δt的初、末代号，初为1，末为2；二维数组槽蓄增量ΔW括号内的第一项为河段上游断面编号，第二项为河段下游断面编号。

因此，按水量平衡方程式计算的全河段累计槽蓄增量∑ΔW水量平衡为：

∑ΔW水量平衡＝Δt(Q(1，1)+Q(1，2)-Q(n，1)-Q(n，2))/2--------------(8)

按动力方程式计算的全河段累计槽蓄增量∑ΔW动力为：

附图说明

图1为本发明的用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法假定糙率系数为0.031和0.035时的1983年洪水杜家台闸计算与实测水位过程线对比图；

图2为本发明的用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法根据水量平衡方程式绘制的图；

图3为本发明的用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法根据动力方程式(2)和(3)绘制的图。

具体实施方式

以下对本发明的用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法的具体实施方式结合图1作进一步详细的说明，但本发明不局限于以下实施例。

实施例一：

根据图1所述的用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法得出的1983年洪水杜家台闸计算与实测水位过程线对比图，分析杜家台闸至周帮河段1983年洪水时的糙率系数。

杜家台闸至周帮河段即洪道长约21km，起止地点分别为杜家台闸和周帮，洪道是引导汉江洪水进入分蓄洪区的通道。

采用1998年实测洪道横断面28个，测量断面编号为1～28，测量精度为IV等水准，断面高程为黄海基面。

杜家台闸下和周帮的观测水位为冻结吴淞基面，换算关系为：冻结吴淞高程＝黄海高程+1.85m。洪水资料时间为1983年10月的8日7时至15日9时，水位观测时段为每小时1次，共计171组数据。杜家台闸下断面的流量过程根据观测的闸上、下游水位按流量公式计算而得。

分析计算杜家台分蓄洪区洪道糙率系数的技术路线是：

首先应假定河段的糙率系数，如在表一中假定的糙率系数为0.035；

按照河段累计槽蓄增量∑ΔW水量平衡和∑ΔW动力应基本相等的原则假定洪道下游末端周帮站的流量，如在表一中1983年10月13日11时假定断面编号为28周帮的流量为2238.14m³/s时，迭代计算累计河段槽蓄增量∑ΔW水量平衡/Δt和∑ΔW动力/Δt均为-351.39(m³/s)，河段27至28的槽蓄增量ΔW动力/Δt为-0.49m³/s。

断面编号小于28大于1的计算断面的流量按照式(7)迭代计算，如在表一中1983年10月13日11时按式(7)计算断面编号27的流量为2238.32m³/s；

断面编号小于28的计算断面的水位按照式(3)迭代计算完成，在表一中1983年10月13日11时按式(3)计算断面编号27的水位为28.12m；

完成上述计算后，若在假定河段糙率系数的情况下，计算的杜家台闸下水位过程线与实测的水位过程线对比误差较小，则可认为该糙率系数即为所求；否则，应重新假定河段糙率系数。

本实施例假定河段糙率系数分别为0.031和0.035，从图1可以看到，当杜家台闸下水位较高，糙率系数为0.031时计算的杜家台闸下水位过程线与实测的水位过程线对比误差较小，糙率系数为0.035时的对比误差较大；当杜家台闸下水位较低，糙率系数为0.031时计算的杜家台闸下水位过程线与实测的水位过程线对比误差较大，糙率系数为0.035时的对比误差较小。洪道糙率系数随杜家台闸下水位涨高而减少的特点符合一般的规律，可以采用。

由于洪道的计算断面有28个，不便全部录入，因此，在表一中只列入了4个计算断面即编号为1、26、27、28，表一中反映的是假定糙率系数为0.035情况下的部分计算结果，节录计算的时间段为1983年10月13日的10时至21时，假定糙率系数为0.031情况下的计算结果也不便录入表一中。

表一：假定n＝0.035时的洪水演进计算结果(节录)

Claims (1)

1.用槽蓄增量关系求解圣维方程组模拟洪水演进的方法，其特征是：将圣维南基本微分方程组在Δt时段简化为水量平衡方程式和动力方程式后，根据水量平衡方程式(1)可计算河段槽蓄增量，根据时段初和末的动力方程式(2)和动力方程式(3)也可计算河段槽蓄增量；

方程组如下：

${\mathrm{\ΔZ}}_1=\frac{{\stackrel{\‾}{Q_1}}^2}{{\stackrel{\‾}{K_1}}^2}L---\left(2\right)$

${\mathrm{\ΔZ}}_2=\frac{{\stackrel{\‾}{Q_2}}^2}{{\stackrel{\‾}{K_2}}^2}L---\left(3\right)$

以上各方程式中：L为河段长度，单位m；Δt为时段，单位s；Z为水位，单位m；Q为流量，单位m³/s；W为河段槽蓄水量，单位m³；ΔW＝(W2-W1)；ΔZ为河段上、下游计算断面的水位差，单位m；

为河段平均流量；为河段平均流量模数；各量脚注1、2分别为时段初、末的时刻；

当已知条件为Q_上1、Z_上1、Q_下1、Z_下1、Q_上2以及Q_下2或已知条件为Q_上1、Z_{上 1}、Q_下1、Z_下1、Q_上2以及Z_下2时，可依据槽蓄增量ΔW的关系求解Z_上2和Z_下2或求解Z_上2和Q_下2：

根据方程式(1)，如假定Q_下2或Z_下2值，可以计算出相应的河段槽蓄增量ΔW；

根据方程式(2)、(3)，继续利用上述假定的Q_下2或Z_下2值，也可以计算出河段槽蓄增量ΔW；

通过计算机程序按上述关系调节Q_下2或Z_下2值，使两种方式计算的河段槽蓄增量ΔW基本相等，此时的Z_上2和Z_下2或Z_上2和Q_下2即为方程组的解；

当计算断面个数为n时，则需要增加n-2个断面流量计算关系式求解方程组；

如以上游起始断面的编号为i＝1，依次往下游的断面编号则为i＝2、i＝3…i＝n-2、i＝n-1、i＝n；

当断面编号1＜i＜n时，根据水量平衡方程式(1)，如断面的流量和河段槽蓄增量用二维数组表示，则断面i-1与断面i之间和断面i与断面i+1之间的水量平衡关系式分别为：

Q(i-1，1)+Q(i-1，2)-Q(i，1)-Q(i，2)＝2ΔW(i-1，i)/Δt---------------(4)

Q(i，1)+Q(i，2)-Q(i+1，1)-Q(i+1，2)＝2ΔW(i，i+1)/Δt---------------(5)

联解方程式(4)和(5)：

Q(i，2)＝(ΔW(i，i+1)-ΔW(i-1，i))/Δt+(Q(i-1，2)+Q(i+1，2))/2+(Q(i-1，1)+Q(i+1，1)-2Q(i，1))/2------------------------------(6)

也可近似写成式(7)等形式：

Q(i，2)＝(ΔW(i，i+1)-ΔW(i-1，i))/2/Δt+(Q(i-1，2)+Q(i+1，2))/2------(7)

式中：二维数组流量Q括号内的第一项为断面编号，第二项为时段Δt的初、末代号，初为1，末为2；二维数组槽蓄增量ΔW括号内的第一项为河段上游断面编号，第二项为河段下游断面编号；

因此，按水量平衡方程式计算的全河段累计槽蓄增量∑ΔW为：

∑ΔW＝Δt(Q(1，1)+Q(1，2)-Q(n，1)-Q(n，2))/2-------------------(8)

按动力方程式计算的全河段累计槽蓄增量∑ΔW为：

$\mathrm{\Σ\ΔW}=\underset{i=2}{\overset{n}{\mathrm{\Σ}}}\mathrm{\ΔW}(i-1,i)---\left(9\right).$

Images