强冲积河流过程二维水沙耦合数值模拟

1、强冲积河流过程二维水沙耦合数值模拟 中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第8期: 962 972 中国科学杂志社 SCIENCE IN CHINA PRESS 岳志远, 曹志先*, 李新, 车涛 武汉大学水资源与水电工程科学国家重点实验室, 武汉 430072; 中国科学院寒区旱区环境与工程研究所, 兰州 730000 * 联系人, E-mail: zxcao 收稿日期: 2007-06-11; 接受日期: 2008-04-18 国家重点基础研究发展计划(编号: 2007CB714106)、国家自然科学基金(批准号: 50459001)和中国科学院知识创新工 程重要方

2、向项目(编号: KZCX3-SW-357-02)资助 摘要 强冲积河流过程泥沙运动非常活跃、河床变形快, 与水流之间存在 强烈的相互作用. 传统的基于简化控制方程的非耦合数学模型违背了基本 守恒律, 只能近似地适用于弱冲积河流过程. 建立普遍适用于强、弱冲积过 程的二维水沙耦合数学模型, 将现有对不可冲刷床面浅水二维流动的、可以 捕捉激波和接触性间断的WAF TVD二阶数值方法扩展至可冲刷床面浅水 二维水沙运动问题. 应用该耦合模型研究了典型冰湖溃决洪水过程. 关键词 冲积河流 洪水 耦合数学模型 泥沙运动 冰湖溃决洪水 近几十年发展了大量的冲积河流数学模型并被广泛应用于研究河流工程、环境、生

3、态与 灾害问题. 但是, 现有模型主要建立在水沙非耦合理论基础之上, 只能近似地适用于输沙强度小、河床变形很慢的弱冲积过程. 然而, 强冲积过程在自然界广泛存在, 其水流急剧变化, 往往诱发非常活跃的泥沙运动和快速河床变形. 高含沙洪水经常发生在中国的黄河以及孟加拉国和巴拉圭等国的一些多沙河流中, 其河床变形速率与水深变化率可能为同一数量级, 是一类典型的强冲积过程. 冲积河流大坝溃决(或拆除)洪水能量大, 必然诱发非常活跃的泥沙运动与显著河床变形1), 也是典型的强冲积过程(如1975年8月特大暴雨导致中国河南省板桥、石漫滩等水库大坝溃决). 冰湖溃决洪水(GLOF, glacier lak

4、e outburst flood)通常发生在高原地区陡峭坡面上, 水流强度大、侵蚀能力强, 可能急剧冲刷坡面, 诱发泥石流灾害, 伴随着全球气候变暖, 许多冰湖具有潜在的溃决危险7, 构成一类典型的强冲积过程. 强冲积过程水流、泥沙与床面之间存在强烈的相互作用, 传统的非耦合数学模拟理论忽略 了其基本力学特征, 其简化的控制方程违背了流体力学守恒定律, 无法适用于强冲积过程. 因此, 需要建立基于完整守恒律和高性能数值格式的水沙耦合数学模拟理论. 国内有学者建立 1) Zech Y, Spinewine B. Dam-break induced floods and sediment move

5、ment-state of the art and need for research. In: First Workshop of EU Project IMPACT. HR Wallingford, 2002 962 中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第8期 非耦合数学模型以研究低水头堤岸溃决水沙过程, 但实质上这仍然局限于弱冲积过程. 近年有学者引入了原本为空气动力学问题而发展的激波捕捉数值方法用于研究不可冲刷床面(定床)溃坝洪水演进过程. Cao等人建立了可冲刷床面(动床)一维水沙耦合数值模型, 清晰地描述了溃坝水流、泥沙及河床变化过程以及相互作用关系, 但

6、局限于一维矩形断面情况. Simpson和Castelltor将Cao等人一维模型扩展至二维, 但其一阶数值格式对于强冲积过程是粗糙的. 本文建立二维水沙耦合数学模型, 将现有应用于定床浅水二维流动的、可以捕捉激波和接触性间断的WAF TVD二阶数值方法扩展至动床浅水二维水沙运动问题, 运用非界面追踪的方法处理干湿边界. 应用该耦合模型研究了典型冰湖溃决洪水过程. 1 数学模型 1.1 控制方程 二维浅水水沙耦合数学模型的基本控制方程包括完整的浑水质量守恒方程和动量守恒方程、泥沙连续方程和河床变形方程, 由流体力学基本守恒律推导. 不失一般性, 这里忽略二阶紊动扩散项. 类似于一维模型12,

7、控制方程可以整理成如下守恒形式: ?U?F?G+=S, (1) ?t?x?y ?h?hu? U=?, (2a) ?hv?hc? hu?2?2hugh/2+?, (2b) F=?huv?huc? hv?huv?, (2c) G=22?hv+gh/2?hvc? (E?D)/(1?p)?2()()()?ghEDu?c?gh(S?S)?sw?0?bxx?2(1?p)?x?, (2d) S=?2?(s?w)gh?c(?0)(E?D)v?gh(Sby?Sy)?2(1yp?)?E?D? ?zD?E=, (3) ?t1?p 963 其中U为守恒向量; F, G为通量向量; S为源项向量; t为时间; x, y

8、为平面坐标; h为水深; u, v为x和y方向的深度平均流速; z为河床高程; c为体积含沙量; g为重力加速度; Sfx, Sfy为x和y方向的阻力坡度; p为床面泥沙孔隙率; E和D为水流底部与河床交界面的泥沙上扬通量和泥沙沉降通量; w和 s为清水和泥沙的密度, 分别取1.0103和2.65103 kg/m3; = w(1?c)+sc为水沙混合体密度; 0=wp+s(1?p)为床沙饱和湿密度; Sbx, Sby为x和y方向的河床底坡. 1.2 封闭模式 应用Manning糙率n计算阻力坡度: h4泥沙沉降通量按下述公式计算: Sx= n2 Sy=n2h43 (4) D=w(1?ca)mc

9、a, (5) 这里w为单颗粒泥沙在清水中的沉降速度; ca为近床体积含沙量, 可以根据平均含沙量计算, 即ca=c, =min(2,(1?p)/c); 指数m=4.45R?0.1, d为泥沙颗粒直径 , R/, 为清水运动黏性系数, 本文取1.110?6 m2/s, s=s/w?1. 对于床面剪切应力较小、水深较大及河床底坡较缓的弱冲积河流, 已经建立了许多泥沙上扬通量的经验公式. 本文研究GLOF事件, 床面剪切应力很大、河床底坡陡峭, 并且在处理干湿边界问题中可能出现水深很小的情况. 从数值计算的稳定性考虑, 选择Zyserman和Fredsoe15 (ZF)公式以估计泥沙上扬通量: ?w

10、(1?ce)mce,c,?E=? (6a) ,?c?0, ce= 0.15226(?c)1.75 0.46+0.331(?c)1.75, (6b) 2/sgd= Shields 参数; c为临界起动Shields参数; u*为床面剪切流速. 这里=u* 2 数值方法 这里将现有应用于定床浅水二维流动的WAF TVD二阶数值方法推广到动床二维水沙耦合问题, 简述如下. 2.1 算子分裂法 应用算子分裂法16离散方程(1): kUip,j=Ui,j?tt(Fi+1/2,j?Fi?1/2,j)?(Gi,j+1/2?Gi,j?1/2), (7a) xy1ppUik,+j=Ui,j+tS(Ui,j),

11、(7b) 这里t为时间步长; x和y为空间步长; i和j为空间节点号; k为时间节点号; p为预报时964 中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第8期 刻节点号; Fi+1/2,j, Gi+1/2,j为x和y方向的数值通量. 河床变形方程(3)离散为 1zik,+j=zik,j+t(D?E)ip,j 1?p. (8) 2.2 数值通量 如下介绍计算方程(7)中x方向界面数值通量16,17Fi+1/2,j计算过程(y方向数值通量的计算与之类似). GGFi+1/2,j=Lwafx,t/2(Ui,j,Ui+1,j), godkUiG,j=Ly,t/2(Ui,j), god

12、kUiG+1,j=Ly,t/2(Ui+1,j), (9a) Fi+1/2,j11N)=(Fi,j+Fi+1,j)?sign(cK)AKFi(+K1/2,j, (9b) 22K=1)(K+1)(K)Fi(+K1/2,j=Fi+1/2,j?Fi+1/2,j, (9c) ?1, rK0,?AK=?(1?cK)2rK (9d) 1,r0,?K?1+rK? ?qi?qi?1?q?q, cK0,?i+1i (9e) rK=?q?q+2+1ii?, cK0,?q?q?i+1i waf这里Lgod y,t/2表示应用Godunov一阶格式计算y方向的中间守恒量, 时间步长取t/2; Lx,t/2表示应用WAF

13、 TVD方法确定x方向的守恒量, 时间步长为t; N为通过该界面上的守恒区间数 )(K)目; Fi(+K1/2,j为第K个波两侧的数值通量差; Fi+1/2,j为黎曼算子. AK为限制函数; 对于q值 本文分别采用如下两种方案: v, c (y方向取u, c, 该情况表示为RV)和h, c (RH)作为q, 并比较了其对数值结果的影响. 在WAF TVD方法中, 应用 HLLC 近似黎曼算子来计算数值通量, 即 ? FL, 0SL,? F, S0S,F*L=FL+SL(U*L?UL),?*L*LHLLCFi+1/2,j=? (10) F , S0S,F=F+S(U?U),*R*RRR*RR?*

14、R ?FR, 0SR. =F(UL,R)为单元左侧和右侧的数值通量; U*L,R为中间守恒量. U*L,R 这里FL,R ?SL,R?uL,R=hL,R?S?S*?L,R?T?c1,Sv?L,RL,R?, (11) *? 965 其中SL, SR, S* 为左、右和中间波波速, aL,R SL=uL?aLqL, SR=uR?aRqR, S*=SLhR(uR?SR)?SRhL(uL?SL). (12) hR(uR?SR)?hL(uL?SL) 类似于Sleigh等人18和Hubbard等人19, 本文采用一种非界面追踪的方法18,19处理干湿边界. 通过上述处理, 模型在时间和空间上具有二阶精度,

15、 其CFL稳定性条件要求柯朗数Cr满足 Cr=maxSx(maxt/x,Symaxt/y1, (13) ) 其中Sx, Sy分别为x和y方向波速. 3 模型性能验证 应用上述耦合数学模型分别进行了定床、无阻力渠道溃坝洪水、复杂定床溃坝洪水和动床溃坝水沙过程的数值模拟试验, 与解析解或现有数值计算结果比较显示本文耦合模型具有较好的精度和适应复杂干湿边界的性能. 作为示例, 这里给出对动床溃坝水沙过程的计算结果。设河道全长150 km, 宽150 m, 坝体位于中间位置(x=75 km), 溃坝前上游水深40 m, 下游水深2 m, 且初始时刻各水力要素沿横向不变. 事实上, 对于该算例而言, 本

16、文二维模型数值结果的各水力要素沿纵向的变化应该与Cao等人一维模型类似, 可以与之比较. 数值计算参数为Cr = 0.95, x = 30 m, y = 30 m, n = 0.03, d = 4.010?3 m, c= 0.045. 对应于ZF泥沙上扬通量计算(6)式, 图1表示本文二维模型与Cao等人12一维模型数值计算结果比较. 从图1看 图1 本文二维模型与Cao等人13 一维模型数值计算结果比较 (a) 水位和河床高程; (b) 流速 966 中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第8期 出, 二维模型与一维模型数值结果符合得很好, 说明本文的二维模型有较好的

17、精度和可靠性, 可以应用于研究强冲积过程中水、沙及河床变形规律. 4 冰湖溃决洪水数值模拟 应用本文二维水沙耦合数学模型研究一冰湖坝体瞬时GLOF演化过程. 分别考虑定床(Cr= 0.01, n=0.05)和动床情况(Cr=0.008, n=0.05, d=10.010?3 m, c=0.06). 数值试验表明当柯朗数Cr?1时计算稳定; 否则, 即使Cr满足对流算子线性稳定性条件(13), 也可能是不稳定的. 这实质上取决于非齐次双曲型方程(1)中源项幅值大小, GLOF所流经区域坡面陡峭、地形异常复杂, 目前尚无普适方法判别其对数值格式稳定性的影响. 图24分别为动床情况下不同时刻GLOF

18、水面高程、含沙量和冲淤厚度分布. 图3显示: GLOF波峰含沙量最大, 而在冰湖以内含沙量很小. 物理上, 波峰水流速度快、床面剪切应力很大, 对床面的侵蚀能力强, 导致含沙量高. 相反, 在湖内水流速度很慢, 以至于无法从湖底起动大量的泥沙, 故含沙量较小. 由图4可以看出, 最大冲刷深度发生在坝体附近, 并且随着时间的推移冲刷严重区域由坝体附近逐渐向湖内发展; 而GLOF波峰冲刷深度较小, 这是因为从上游输送下来的大量泥沙会逐渐沉积到河床上. 一般而言, 含沙量与床面冲刷深度沿洪水主流方向的分布与Cao等人一维模型所揭示的规律是一致的. 动床与定床条件下GLOF演化过程是有显著差异的, 导

19、致GLOF波峰到达的时间及水位变化过程不尽相同. 图5显示冰坝下游L=4.1 km和L=8.9 km处动床与定床条件下水位变化过程. 由图5可以看出, 动床情况下GLOF波峰到达时间均滞后于定床. 与定床情况相比, 动床GLOF对河床的冲刷使其损耗了部分能量, 并且河床形态变化必然引起形态阻力增加, 这两方面的综合效应促使GLOF传播速度减小; 另一方面, 动床GLOF挟沙水流密度增加, 重力作用增大, 从而成为洪水传播的驱动力. 与ZF泥沙上扬通量计算(6)式对应, 动床GLOF冲刷河床能量损耗与形态阻力增加之综合效应超越了挟沙水流密度增加而引起的重力作用变化之效应, 致使动床GLOF传播速

20、度较慢, 波峰到达时间滞后. 从图5中还可以看出, 在动床与定床条件下GLOF最大水位也有显著差异, 床面冲刷对水位变化过程的影响是显著的. 图6为动床条件耦合解、动床条件非耦合解和定床情况下GLOF波峰与坝址之间距离随时间的变化. 在3种条件下, 冰坝溃决后在1500 s时刻GLOF分别传播了约13.3, 15.7及17.4 km. 这首先显示出定、动床条件下GLOF传播速度的差异. 不仅如此, 在定、动床条件下GLOF传播速度远远快于一般平原河道洪水传播速度, 这是因为GLOF所流经的高原区域坡面陡峭、重力作用相对于床面阻力而言占主导地位. 图6还显示出动床条件下耦合解与非耦合解之间存 在

21、显著差异. 这说明在GLOF传播过程中, 水流、泥沙及河床变形三者之间存在着强烈的相互作用. 图7表示动床与定床情况下坝址附近流量过程线. 可以看出, 坝址处动、定床条件下流量过程差异显著. 一般而言, 动床流量峰值大于定床情况, 这主要基于两方面原因: 一方面, 坝址处溃决洪水能量最大, 坡面所遭遇的冲刷最为严重, 一定程度上增大了坝址附近的坡面比降; 另一方面, 动床情况的挟沙水流重度大于定床情况的清水水流. 两方面因素的作用使得动床条件GLOF重力变化之效应超越了能量损耗和形态阻力增加之效应, 使得动床情况在坝址 967 图2 动床情况下在不同时刻GLOF水面高程分布 (a) t = 0

22、 s; (b) t = 300 s; (c) t = 1500 s 968 中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第8期 图3 动床情况下在不同时刻的GLOF含沙量分布 (a) t = 300 s; (b) t = 1500 s 图4 动床情况下在不同时刻的GLOF冲淤厚度分布 (a) t = 300 s; (b) t = 1500 s 处具有更大的峰值流量. 床面泥沙上扬的定量化是水沙耦合数学模拟理论所涉及的关键问题之一. 然而, 目前对强冲积过程泥沙运动机理的认识远未清楚, 尚无普遍适用的强冲积过程泥沙上扬通量计算模式. 当泥沙上扬通量增加时, 如对应于Cao等人2

23、0所应用的泥沙上扬通量公式, 其动床GLOF传播速度则明显快于定床情况. 这与本文应用方程(6)所揭示的规律(图5和6)不同, 这是由于动床GLOF冲刷河床能量损耗与形态阻力增加之综合效应尚不足以超越挟沙水流密度增加而引起的重力作用变化之效应. 毫无疑问, 深入研究强冲积过程水流、泥沙运动与河床相互作用机理, 建立适用于强冲积过程的泥沙上扬定量化模式是十分必要的. 目前的实际应用中, 建议用实际水文观测资料率定泥沙上扬定量化模式参数和阻力参数等, 再将本文建立的水沙耦合模型应用于相同或相似环境的冲积河流过程. 图8为分别采用RV和RH的模型数值结果比较. 从图8中可以看出, 在计算限制函数时,

24、 由于q值选取方法的不同, 模型数值结果(水位、河床、流速、含沙量变化过程)之间的差异很 969 图5 动床与定床情况下GLOF水位变化过程 (a) L = 4.1 km; (b) L = 8.9 km 图6 动床与定床情况下GLOF波峰与冰坝的距离 大, 这在一定程度上也表明数值格式的差异能够影响数值计算结果. 显然, 为了进一步检验模型的精度和可靠性, 迫切需要实测资料对模型进行验证, 以将其推广应用于实际. 970 中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第8期 图7 动床与定床情况下GLOF在坝址处流量过程线 图8 分别采用RV和RH的模型数值结果比较 (a) 水

25、位及河床高程过程; (b) x方向流速过程; (c) y方向流速过程; (d) 含沙量过程 971 5 结论 建立了可冲积床面浅水二维水沙耦合数学模型, 将现有对不可冲积床面浅水二维流动的WAF TVD二阶有限体积数值方法推广用于数值求解可冲积床面浅水二维水沙耦合数学模型的双曲型控制方程, 普遍适用于强、弱冲积河流过程. 运用该模型研究了冰坝瞬时溃决洪水传播过程及其诱发的泥沙运动和河床演化. 由于GLOF流经区域坡面坡度大, 重力作用占主导地位, 导致GLOF传播速度很快. 可冲积床面上GLOF冲刷河床能量损耗、形态阻力变化以及挟沙水流密度增加而引起的重力作用变化等效应之综合作用, 可能使得G

26、LOF传播速度显著地不同于类似条件下不可冲积床面GLOF的传播, 显示强冲积河流过程河床变形与水流之间强烈的相互作用. 特别地, 必须发展耦合数学模拟理论以尽可能充分地反映水流、泥沙运动和河床演化过程. 深入研究强冲积过程泥沙上扬机理、建立合理的定量化模式是十分必要的. 另外, 数值格式的差异也会影响模型计算结果, 这就需要以实测资料对模型进行进一步的检验. 参考文献 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Zhang R J, Xie J H. Sedimentation Research in China-Systematic

27、Selections. Beijing: China Water and Power Press, 1993 Cenderelli D A, Wohl E E. Flow hydraulics and geomorphic effects of glacial-lake outburst floods in the mount Everest region, Nepal. Earth Surf Proc Landf, 2003, 28: 385407DOI Cenderelli D A, Wohl E E. Peak discharge estimates of glacial-lake ou

28、tburst floods and “normal” climatic floods in the mount Everest region, Nepal. Geomorphology, 2001, 40: 5790 OConnor J E, Baker V R. Magnitudes and implications of peak discharges from glacial lake Missoula. Geol Soc Am Bull, 1992, 104: 267279DOI Baker V R, Benito G, Rudoy A N. Paleohydrology of lat

29、e Pleistocene superflooding, Altai Mountains, Siberia. Science, 1993, 259(5093): 348350 Liu S Z, Li H X, Yan Y, et al. Assessment of bursting hazards of the ice lakes in Luozha county, Tibet. J Mt Sci, 2003, 12(1): 128132 Che T, Li X, Mool P K, et al. Monitoring glaciers and associated glacial lakes

30、 on the east slopes of mount Xixia-bangma from remote sensing images. J Glaciol Geocryol, 2005, 12(6): 801805 梁林, 倪晋仁. 黄河溃堤过程数学模型及其模拟方法. 中国科学E辑: 技术科学, 2002, 32(5): 618627 Zhang X Z, Wang G Q. Flow analysis and scour computation of the dike-break. J Sed Res, 2002, 1: 1824 Zoppou C, Roberts S. Numeric

31、al solution of the two-dimensional unsteady dam break. Appl Math Model, 2000, 24: 457475 Valiani1 A, Caleffi V, Zanni A. Case study: Malpasset dam-break simulation using a two-dimensional finite vol-ume method. J Hydraul Eng, 2002, 128: 460472DOI Cao Z X, Pender G, Wallis S, et al. Computational dam-break hydraulics over erodible sedim