《强冲积河流过程二维水沙耦合数值模拟.doc》由会员分享，可在线阅读，更多相关《强冲积河流过程二维水沙耦合数值模拟.doc（16页珍藏版）》请在金锄头文库上搜索。
1、强冲积河流过程二维水沙耦合数值模拟 中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第8期: 962 972 中国科学杂志社 SCIENCE IN CHINA PRESS 岳志远, 曹志先*, 李新, 车涛 武汉大学水资源与水电工程科学国家重点实验室, 武汉 430072; 中国科学院寒区旱区环境与工程研究所, 兰州 730000 * 联系人, E-mail: 收稿日期: 2007-06-11; 接受日期: 2008-04-18 国家重点基础研究发展计划(编号: 2007CB714106)、国家自然科学基金(批准号: 50459001)和中国科学院知识创新工 程重要方向项目
2、(编号: KZCX3-SW-357-02)资助 摘要 强冲积河流过程泥沙运动非常活跃、河床变形快, 与水流之间存在 强烈的相互作用. 传统的基于简化控制方程的非耦合数学模型违背了基本 守恒律, 只能近似地适用于弱冲积河流过程. 建立普遍适用于强、弱冲积过 程的二维水沙耦合数学模型, 将现有对不可冲刷床面浅水二维流动的、可以 捕捉激波和接触性间断的WAF TVD二阶数值方法扩展至可冲刷床面浅水 二维水沙运动问题. 应用该耦合模型研究了典型冰湖溃决洪水过程. 关键词 冲积河流 洪水 耦合数学模型 泥沙运动 冰湖溃决洪水 近几十年发展了大量的冲积河流数学模型并被广泛应用于研究河流工程、环境、生态与
3、灾害问题. 但是, 现有模型主要建立在水沙非耦合理论基础之上, 只能近似地适用于输沙强度小、河床变形很慢的弱冲积过程. 然而, 强冲积过程在自然界广泛存在, 其水流急剧变化, 往往诱发非常活跃的泥沙运动和快速河床变形. 高含沙洪水经常发生在中国的黄河以及孟加拉国和巴拉圭等国的一些多沙河流中, 其河床变形速率与水深变化率可能为同一数量级, 是一类典型的强冲积过程. 冲积河流大坝溃决(或拆除)洪水能量大, 必然诱发非常活跃的泥沙运动与显著河床变形1), 也是典型的强冲积过程(如1975年8月特大暴雨导致中国河南省板桥、石漫滩等水库大坝溃决). 冰湖溃决洪水(GLOF, glacier lake o
4、utburst flood)通常发生在高原地区陡峭坡面上, 水流强度大、侵蚀能力强, 可能急剧冲刷坡面, 诱发泥石流灾害, 伴随着全球气候变暖, 许多冰湖具有潜在的溃决危险7, 构成一类典型的强冲积过程. 强冲积过程水流、泥沙与床面之间存在强烈的相互作用, 传统的非耦合数学模拟理论忽略 了其基本力学特征, 其简化的控制方程违背了流体力学守恒定律, 无法适用于强冲积过程. 因此, 需要建立基于完整守恒律和高性能数值格式的水沙耦合数学模拟理论. 国内有学者建立 1) Zech Y, Spinewine B. Dam-break induced floods and sediment movemen
5、t-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)mca,
9、(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(1?
10、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), go
12、dkUiG+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表示应用WA
13、F 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处理干湿边界. 通过上述处理, 模型在时间和空间上具有二阶精度, 其CFL稳定性条件要求柯朗数Cr满足 Cr=maxSx(maxt/x,Symaxt/y1, (13) ) 其中Sx, Sy分别为x和y方向波速. 3 模型性能验证 应用上述耦合数学模型分别进行了定床、无阻力渠道溃坝洪水、复杂定床溃坝洪
