用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法 技术领域 本发明涉及防洪减灾科学技术领域， 特别是根据水量平衡方程式与动力方程式的 槽蓄增量关系求解圣维南方程组模拟洪水演进的方法。
背景技术 我国是一个洪涝灾害频繁发生的国家， 全国 1/2 的人口， 1/3 的耕地和主要大城市 都处于江河洪水位以下， 受洪水威胁的地区工农业产值占全国 2/3。 为了达到减轻洪水灾害 以及利用水资源为国民经济与社会发展服务的目的， 建立模拟洪水演进的数学模型， 在河 道、 渠道、 分蓄洪区和水库合理确定防洪工程规模、 洪水预报、 科学调度洪水等方面发挥着 重要作用。
描述河渠非恒定流运动规律的一阶拟线性双曲型偏微分方程， 在 1871 年首先由 法国学者圣维南建立。在一般情况下， 无法求出其普遍的解析解。目前， 只能使用近似的计 算方法求解。近似的计算方法大致可分为两类 ：
1、 水力学方法 ：
直接差分法就是其中的一种方法。直接差分法是用偏差商代替偏微商， 将基本微 分方程离散化为差分方程， 求在自变量域 x ～ t 平面差分网格上各节点近似数值解的方法。 基本微分方程中的偏微商可用不同形式的偏差商即差分格式代替， 从而得到不同的差分方 程。 但从目前的技术手段而言， 研究成果一般用于规划设计阶段， 在流域水文预报模型中尚 未采用。
2、 水文学方法 ：
水文学方法如著名的 “马斯京根法” ， 是将基本微分方程组分别简化为水量平衡方 程式和近似为蓄泄方程式， 用差分形式合解水量平衡方程式与蓄泄方程式。 缺点是 ： ①概化 相对较粗， 难于详细分析河道内复杂水流情势以及主要特征 ； ②难于提供河段内水位的时 空分布和动态模拟洪水在河段内的演进过程。
发明内容 本发明的用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法， 是将圣维南 基本微分方程组在 Δt 时段简化为水量平衡方程式和动力方程式后， 根据水量平衡方程式 和动力方程式可以分别计算河段槽蓄增量 ： 一是根据水量平衡方程式可计算槽蓄增量， 二 是根据时段初、 末的两个动力方程式也可计算槽蓄增量。由于两种方法计算的河段槽蓄增 量应相等， 因此， 这一发明就为水力学方法求解圣维南方程组模拟洪水演进提供了一种新 的、 可靠的方法。
本发明的用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法， 于 2009 年 11 月下旬已在湖北省水利水电勘测设计院进行的杜家台分蓄洪区可行性研究报告修编工作 中得到了初步应用， 通过模拟洪水演进计算， 结果发现没有必要兴建原报告的支洪道项目 ( 含出江口控制闸 )， 若实施， 可节省投资 2 亿余元。
本发明的用槽蓄增量关系求解圣维南方程组模拟洪水演进的方法， 属于水力学方 法范畴。本方法既与水文学方法和现行的水力学方法有相似的地方， 也有不同的地方。
与水文学方法相似的地方是都涉及到了河段槽蓄内容 ； 不同的是水文学方法如 “马斯京根法” 是将动力方程式近似变为蓄泄方程式后与水量平衡方程式联解， 而本方法是 直接将水量平衡方程式与动力方程式联解。
与现行的水力学方法相比， 相同之处是都涉及到了需要河道地形资料和依据水文 资料率定河段糙率 ； 不同的是， 本方法将基本微分方程在 Δt 时段简化为水量平衡方程式 和时段初、 末的两个动力方程式后联解。
具体方法如下 ：
现将圣维南基本微分方程组在 Δt 时段简化为水量平衡方程式 (1)、 动力方程式 (2) 和动力方程式 (3) ：
以上各式中 ： L 为河段度长 (m) ； Δt 为时段 (s) ； Z 为水位 (m) ； Q 为流量 (m3/s) ； q为区间流量 (m3/s) ； W 为河段槽蓄水量 (m3) ； Δ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( 动力 ) 是在时刻 t1 的水面线 Z 上 1 ～ Z 下 1 和时刻 t2 的水面线 Z 上 2 ～ Z 下 2 所包围的水体， 即由上、 下游计算断面 分别在时刻 t1 和 t2 的水位转换成过水断面积的差值平均后乘以河段长度得出的体积。
当已知条件为 Q 上 1、 Z 上 1、 Q 下 1、 Z 下 1、 Q 上 2 以及 Q 下 2( 或 Z 下 2) 时， 可依据槽蓄增量 ΔW 的关系求解 Z 上 2 和 Q 下 2( 或 Z 下 2) ：
根据方程式 (1)， 如假定 Q 下 2( 或 Z 下 2) 值， 就可以计算出相应的河段槽蓄增量 ΔW( 水量平衡 )。
根据方程式 (2)、 (3)， 继续利用上述假定的 Q 下 2( 或 Z 下 2) 值， 就可以计算出河段槽
蓄增量 ΔW( 动力 )。
通过计算机程序按上述关系调节 Q 下 2( 或 Z 下 2) 值， 使两种方式计算的河段槽蓄增 量 ΔW 基本相等， 此时的 Z 上 2 和 Q 下 2( 或 Z 下 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 时假定断面编号为 3 28( 周帮 ) 的流量为 2238.14(m /s) 时， 迭代计算累计河段槽蓄增量∑ ΔW( 水量平衡 )/ Δt 和∑ ΔW( 动力 )/Δt 均为 -351.39(m3/s)， 河段 27 至 28 的槽蓄增量 ΔW( 动力 )/Δt 为 -0.49(m3/s)。
断面编号小于 28( 周帮 ) 大于 1( 杜家台闸下 ) 的计算断面的流量按照式 (7) 迭代 计算， 如在表一中 1983 年 10 月 13 日 11 时按式 (7) 计算断面编号 27 的流量为 2238.32(m3/ 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 时的洪水演进计算结果 ( 节录 )