水 利 学 报
SHUILI XUEBAO
第5期
1998年4月
黄河下游河道水沙运动仿真模型的开发研究
程晓陶 黄金池
(中国水利水电科学研究院)
(黄河水利委员会河务局)
薛云鹏
摘 要 本文介绍的黄河下游宽浅河段汛期水沙运动仿真模型建立在二维不恒定流与不平衡输沙理论的基础上,采用无结构不规则网络,结合了有限差分法和有限体积法的优点,能较充分地反映滩区、生产堤与滩槽冲淤变化对洪水演进的显著影响,并具有实时动态图象显示和人机对话功能.关键词 二维不恒定流,水沙运动仿真.
黄河下游滩宽槽浅,汛期来水来沙变幅很大.大堤间滩地面积占河道面积83%,有居民近140余万.滩区生产堤纵横,在纵比降为1/6000的河道内,滩面横比降可达1/2000;平滩流量下降到约4000m3/s.洪水一旦破堤漫滩;顶冲大堤,顺堤行洪,甚至倒灌,对滩区人民生命财产及黄河大堤的威胁更大.漫滩行洪和滩槽水沙激烈交换可引起洪峰严重变形,且河槽冲淤幅度与水位涨落达到同等
量级.大幅快速的河床调整又反过来影响水流的挟沙能力与河道的行洪能力,引起水位-流量关系的复杂变化,使洪水预报异常困难.
对这种大范围的、在计算域内部受堤防溃决等突发性因素影响十分显著、水沙交换激烈的动边界、动床洪水现象,一般一维或一维扩展模型很难模拟.已有的洪水预报模型或河势演变模型对漫滩洪水大都作了简化,难以充分体现生产堤溃决、滩地行洪与滩槽水沙交换对洪水传播变形的影响.新研建的黄河下游河道水沙运动仿真模型,既要求计算精度高,以便为滩区洪水风险分析提供科学依据;又要求运算速度快,在大洪水发生时能及时准确为滩区居民发布洪水预警报信息.如此矛盾的要求增加了模型设计的难度.
目前,黄河下游夹河滩至高村(包括封丘倒灌区)防洪险要河段的水沙运动仿真模型已通过鉴定,其主要特点包括:①以一、二维不恒定流与不平衡输沙理论为基础,采用无结构不规则网格,可以较好地概化河道主槽、滩地的地形特征,充分反映生产堤、渠堤等对水沙运动的影响.②泥沙计算在黄金池方法的基础上作了改进,引入了张红武的一套挟沙能力公式,使黄河下游主槽冲淤变化对洪水演进的影响能够较好体现出来,并能适应高含沙洪水的计算.③模型具有在运算过程中实时动态显示流量、水位、洪水淹没范围或泥沙冲淤计算成果的功能.④模型中设计了人机对话的窗口,对生产堤在洪水期间是否发生决口这种带有随机性的问题,模型使用者有机会根据实际情况或经验进行选择.⑤模型运算速度显著提高,对1427.8km2的计算域模拟6天洪水过程,仅需十余分钟机时,可满足为洪水预报服务的要求.⑥开发了模型的菜单对话式图形信息管理系统,使模型所需大量基本信息的前处理与计算成果的后处理工作大为简便.⑦计算结果与损失评估模型相结合,可以快速算出相应洪水条件下的受灾人口及可能造成的经济损失.
[1]
1 水沙运动的计算方法
1.1 基本控制方程 从黄河水沙运动特点出发,模型中选垂向平均单宽流量Q与水深H作为基本状态变量,其深度积分形式的控制方程组可记作如下的质量与动量守恒形式:
本文于1996年7月1日收到.
—12—H+ M+ N=q,水流连续方程 t x y
M (uM) (vM) Zn2uu2+v2水流动量方程 +++gH+g=0,
t x y xH1/3
(1)(2)
N (uN) (vN) Zn2vu2+v2+++gH+g=0,(3) t x y yH1/3
式中:H—水深;Z—水位;q—源汇项,本模型中代表有效降雨强度,不计降雨情况下q=0;M,N—分别为x、y方向的垂向平均单宽流量;u,v—分别为垂向平均流速在x、y方向的分量;n—曼宁糙率系数;g—重力加速度.
基本状态变量与其他交量之间有以下转换关系:
Z=H+B,u=M/H,
v=N/h,
式(4)中B为地面高程.
在考虑泥沙冲淤计算时,基本控制方程组中要增加泥沙运动方程与河床变形方程.根据不平衡输沙理论,这两个方程通常可表达为如下形式:泥沙运动方程
(HS) (HuS) (HvS) (+ 2(HS)=aω++-es高含(S＊-S),2HS) t x y x y
2
2
(4)(5)
(6)
(7)
B河床变形方程 r′=αω(S＊-S),(8)
t
式中:S,S＊—分别为深度平均含沙量与水流挟沙力;ω—悬移质泥沙深度平均沉速;α—悬移质泥沙恢复饱和系数;es—悬移质泥沙扩散系数;γ′—淤积物干容重;其余符号与水流模型意义相同.1.2 水流控制方程的简化与离散方式 黄河滩区计算域大,河势复杂多变,控导工程、生产堤、渠堤密布其间.在规则网格情况下,若为省计算时间而将网格空间步长放大,则会引起严重的地形失真,难以满足洪水风险评价对计算精度的要求.若加密网格以提高对地形的分辩率,又必然会增大计算量,难以满足缩短运算机时、为洪水预报服务的要求.
为此,本模型在开发中保留了体积积分形式的有限差分法的优点[2],如在网格周边通道上计算流量、泥沙通量,在网格的形心计算水深、含沙量、挟沙力与底床的冲淤变化;同时又吸收了有限体积法中非结构性不规则网格的优点[3].通道上定义的单宽流量,不限于正交的X或Y的方向,而可以是任意的通道法线方向的流量,同时,水深与流量在时间轴上分层布置,交替求解,物理意义清楚,并有利于提高计算稳定性,如图1所示.
a.H与Q的空间布置方式b.H与Q的时间交错计算方式
图1 基本状态变量的离散化布置方式
将式(5),(6)代入式(1),则水流连续方程式改写后对计算域进行面积分,再利用高斯定得,可得
—13—∫A
HdA+(Hu·n)dl= t
∮qdA,∫
A
(9)
式中:u为计算域边缘上任一点的流速矢量,n为该点的外法线方向单位向量.令Q=Hu·n,则Q为
任意n方向垂向平均单宽流量,当n取为x、y方向的单位向量时,Q即为M、N.
由于黄河下游河道比降较小,滩区地势平坦开阔,洪峰流量靠河宽变化调整,水位变幅不大,可以假定水深H随时间的变化在一个有限的网格内差异很小,降雨在一个网格内也是均匀分布的,则式(9)简化成:
HA t+
Qdl=Qq.∮
L
K
(10)
对任一K边形网格,上式左端第2项的线积分可记为:
Qdl=∑QL∮
L
k=1
kk
,(11)
式中,A为网格面积,L为通道长度,下标k为对K边形网格逆时针方向的通道编号.
设对任意网格,Q流入为正,流出为负,则连续方程对i网格的显式离散化形式可以写成
HiT+2DT
式中下标i为网格的编号.
式(12)表明,在已知T时刻水位的情况下,为求出网格中T+2DT时刻的水位,只要能合理确定各通道上T+DT时刻水流的单宽通量Q既可.为此,针对黄河特点,将通道分为河道型、滩地
型、有堤型等,采用不同假设进行简化计算.
河道型通道设在主槽中.对于黄河下游平原型河道,动量方程中仅忽略了对流项,即假定第j通道上的垂向平均单宽流量Qj不仅与相邻两网格的水位Zj1,Zj2有关,而且与前一时刻流量有关,则动量方程的离散化形式简化为:
+DTQTj
=HTi
2DT+
Ai
k=1
+DT
Likik∑QT
K
+2DTqT+DT,
(12)
=QjT-DT
Zj2-Zj1n2QQ-2DT×ghj-2DT×g,/3
DLjh7j
T
T
TT
(13)
式中hj为j通道上的平均水深,一般由下式计算:
hj=DCj1Hj2+DCj2Hj1/DLj,
有(限14)
式中DLj=DCj1+DCj2,为空间步长,如图2(a)所示.
考虑到黄河河道宽浅,有时需将一个横断面分成若干网格,而网格间高差较大,成阶梯状,式(14)计算的hj可能过大.因此,判断河底比降大于0.0006时,hj改由下式插值给出
TT
hj=DCj1Zj2-Bm+DCj2Zj1-cBkm
gpt/DL,ja
g
(15)
式中Bm为相邻网格中较高的河底高程.
滩地型通道设在陆地上.漫滩洪水传播主要受重力与阻力的作用,可按扩散波运动进行处理.从动量方程中完全略去惯性项,其离散化形式简化为:
to
T+DTQj
,
=sign
ZjT1-ZTj25/3
hj
T
ZTj2-Zj1
DLj
1/2
1·n,
(16)
式中DL的定义为相临网格形心距在通道法线方向的投影,如图2(b)所示.
当滩地水深达到一定深度后,顺滩行洪状况型等同于河道主槽.模型计算中将根据给定的判别条件由式(16)自动转换为式(13)的主槽算法.
对有阻水建筑物的通道,可分别概化成连续堤防、有缺口堤防及有桥涵的堤防等几种情况,采用适当的堰流或孔口出流等经验公式进行计算.
根据模拟对象的特点,在网格类型上也区分了主槽、嫩滩、老滩.网格的非结构性设计,使其既便于针对主槽的游荡性进行计算域的调整,又便于今后计算域向花园口至陶城铺河段扩张;网格的不—14—规则性,使地形概化更接近实际,能更好反映复杂的河势变化与工程条件对洪水演进的影响;网格分类上的灵活性,使其能针对计算对象的基本特点采取不同的简化方式,既抓住主要矛盾,又节省计算时间.
a.河道型通道
图2 不规则网格空间步长的分类处理
b.滩地型通道
为适合黄河的特点,模型开发中对水流与泥沙计算的结合、网格形心流速的算法、滩唇的处理方
法、出口断面无水位-流量关系及滩区横比降引起的数值效应等问题均作了有益的探讨,受篇幅限制,详细计算与处理方法将另文介绍,或请参见文献[4].
2 模型的验证
在模型开发阶段,先后以1982年与1973年两场洪水资料对模型进行了验证.前者是三门峡水库建成后黄河下游发生的最大一场洪水,代表了大流量、小含沙的情况;后者则代表了小流量、高含沙的情况.模型经多种验证均得到令人满意的结果[4,5].高含沙洪水的计算及河床冲淤成果的详细分析,笔者将另文介绍,本文仅以1982年洪水验证为例.
模型验证范围取黄河下游夹河滩至高村河段,河道长83km,大堤间距最宽达20km,平均约11km,面积908km2,包括长垣、东明两个大滩,面积均在200km2以上.滩区生产堤完整,两岸生产堤堤距约3km左右.验证结果表明:①洪水淹没范围与实际相当吻合,而计算图更清楚地给出了整个计算域上的最大水深分布.计算成果可对任意时刻绘制出洪水平面二维淹没水深分布图与相应的流速矢量分布图[4];②合理模拟了洪水破堤漫滩后顶冲大堤,顺堤行洪、甚至倒灌的现象.比较图3、图4,前者较好模拟出了洪峰传播变形的特点;而后者由于洪水不能进滩,四个断面的洪水波形几乎不变,后期出口断面水位得不到滩地归槽洪水的补充,明显偏低;③计算最高水位在整个域中分
图3(a) 1982年高村站计算与实测流量过程比较图3(b) 1982年高村站计算与实测水位过程比较
布合理(图5);④滩槽冲淤分布合理,充分体现了大洪水期间“淤滩刷槽”的现象,主槽中出现了一条清晰的连续的冲刷槽,反映了黄河汛期“大水出好河”的特点.根据实测资料,“82.8”洪水期间,泥沙淤积主要分布在生产堤间嫩滩上,约占滩地淤积量的95%,而生产堤外滩地的淤积量很少,计
—15—算结果表明泥沙淤积也是绝大多数分布在生产堤间的滩地上.主槽与嫩滩冲淤幅度的量级是合理的,主槽冲刷深度约为1.0—2.5m,嫩滩淤积厚度约为0.2—1.0m.
图4(a) 假定生产堤不溃的流量过程模拟结果图4(b) 假定生产堤不溃高村站水位过程的比较
图5(a) 1982年洪水黄河北岸大堤计算最高水面线与实测资料比较
图5(b) 1982年洪水黄河南岸大堤计算最高水面线与实测资料比较
3 模型的实时动态显示功能
以往的数值计算模型总是先“盲算”,只有待计算结束后,再进行数据分析,才能知道计算结果.而黄河水沙运行模型计算域广、数据量大,为了改变模型调试、运行效率低的状态,并适应计算中人-机对话的需要,专门开发了一整套绘图子程序,与计算程序一起运行,使各种计算结果能够实时动态显示出来,包括:①以不同颜色表示的水深或泥沙冲淤深度的平面分布.两者之间在计算中可以随时用指定的功能键进行往复切换,想看水深分布还是冲淤分布,随意选择.②以不同颜色线段表示的若干指定河道断面的流量变化过程与若干指定地点的水位变化过程.在有实测资料的情况下,还可以显示计算结果与实测过程的比较。③一些文字信息,如水量平衡系数、计算时间步长、计算时间步数、以年月日时表示的计算时间等.④平面图的比尺.当计算域调整时,只需改变图的中点坐标和比尺系数,即可调整模型的位置和大小,同时显示的平面图的比尺也自动调整.
4 模型的人机对话功能
以往的数学模型,在给定初始条件与边界条件之后,就成为一个确定性的模型,只能一算到底.但黄河的情况非常复杂,尤其生产堤漫溃、冲溃或抢险堵口,具有很大的随机性.为此专门设计了模型的人机对话框,包括自动对话与强制对话两种形式:①自动对话.当堤防出现洪水漫顶时,对话框里自动报警,显示该堤防所在通道的编号(同时平面图上该段通道改变成亮红色,以提示其位置)、堤顶高程、两侧网格的编号、地面平均高程与水位,以及推荐的缺口宽、缺口深等.决策者可以根据—16—经验和当时情况决定堤防溃或不溃,溃的话是全溃还是部分溃,部分溃的话是否接受推荐值,不接受推荐值的话,重新输入缺口宽、缺口深等.②强制对话.实际防汛过程中存在“扒口分洪”、“抢险堵口”、“堤防加高”、“非漫顶的溃决”等多种情况.计算中需要时可以按下强制对话功能键,通过对话框输入所关心的通道号.根据其位置提示,如发现有误,可以重输.确认后,则如同自动对话一样,计算机将提供有关信息,决策者可以改变堤防状态.强制对话亦可仅用于信息查询.
5 结束语
模型经过检验之后,1995年汛前按最新河道大断面及生产堤资料调整了模型,并将计算域扩展到包括封丘倒灌区的1427.8km2,模拟计算了4种不同量级的洪水,即1976年型花园口站实测9210m3/s,1982年型花园口站实测15300m3/s洪水,1958年型花园口站22300m3/s设计洪水以及1982年型花园口站30000m/s设计洪水.研究了在现状河道地形与生产堤废除条件下,历史典型洪水重现或规划最大洪水发生时滩区的淹没情况及夹河滩至高村河段、尤其是封丘倒灌区的滞洪消峰作用,并为进行损失评估和减灾对策研究提供了水情依据.
目前,模型正在向黄河防洪减灾计算机网络系统移植.1997年汛前,模型的计算域将进一步向上下游扩展至花园口～孙口河段.
致谢 模型的开发工作是在与黄委会河务局、河南省河务局的密切合作下完成的,并得到了黄委会吴致尧、赵业安、刘月兰、潘贤娣、张红武等众多专家的关心和指点;陈喜军为模型开发了图像信息处理系统,方便了模型的使用.谨此深表感谢!
3
参考文献
1 张红武,等.黄河高含沙洪水模型相似律的研究.黄委会水科院,1993,12.2 刘树坤,等.全民防洪减灾手册.沈阳:辽宁人民出版社,1993.
3 谭维炎,胡四一.二维浅水流动的一种普适的高性能格式有限体积Osher格式.水科学进展,1991,2(3).4 程晓陶,等.黄河下游滩区水沙运行数学模型及减灾措施的研究.“八五”国家重点科技攻关报告(85-926
-01-02-01),1995,12.
5 程晓陶,杨磊,等.分蓄洪区洪水演进数值模型.自然灾害学报,1996,5(1).
2-Dsimulationmodeloffloodroutingandsedimenttransportinthe
lowerreachoftheYellowRiver
ChengXiaotaoHuangJinchi
(ChinaInstituteofWaterResourcesandHydropowerResearch)
XueYunpeng
(RiverBureauoftheWaterResourcesCommitteeoftheYellowRiver)
Abstract Basedonthetheoryofunsteadyflowandnon-equilibriumsedimenttransport,a2-DsimulationmodeloffloodroutingandsedimenttransportinthelowerreachoftheYellowRiverhasbeendeveloped.Themodelisbuiltbynon-structuralirregulargrids,combiningwithfinitedifferenceandfinitevolumemethods,inwhichthesignificantinfluencesofoverbankstorage,localleveandscouring(silting)intheprocessoffloodpropagationcanbewellconsidered.
Keywords 2-Dunsteadyflow,floodrouting,sedimenttransport,simulation.
—17—
因篇幅问题不能全部显示，请点此查看更多更全内容