堤防溃决的流动及冲刷坑计算分析.docVIP

堤防溃决的流动分析及冲刷坑计算 张修忠，王光谦（清华大学 水沙科学教育部重点实验室） 摘要：以平面二维浅水方程和不平衡输沙理论作为建模基础，采用高分辨率有限元格式捕捉溃坝涌波的传播，建立了堤坝瞬间溃决后水流运动和河床冲刷的二维数学模型。计算了不同溃口宽度、不同流量和不同粒径土体组成条件下的堤坝瞬时溃决的水流演进和冲坑发展，为钢木土石组合坝的结构受力及稳定性分析提供水动力学条件。 关键词：堤防溃决； 冲刷坑； 高分辨率有限元格式； 数值模拟 基金项目：国家自然科学基金及水利部联合资助重大项目，(项目编号。作者简介：张修忠(1972-)，男，山东临沂人，博士生。 1 引言 河道堤防的溃决往往给人民的生命和财产造成巨大的损失。例如，1996年8月河北省滹沱河堤防溃决164m，由决口倾泄而出的2.3亿立方米的洪水淹没了23个村庄和27万亩农田[1]。又如1998年长江九江大堤决口不但加重了洪水的灾害，而且使抗洪工作面临更加严峻的局面。决口的快速封堵成功与否以及封堵费用的高低在很大程度上取决于决口水流的水力特性和溃口附近冲坑的发展变化。因此，对堤防溃决问题的深入研究对于防灾减灾以及灾害的治理有着十分重要的现实意义。 堤防溃决的原因主要有两种[2]，一是管涌发展到一定程度导致堤身塌陷发生漫流破坏，二是洪水位超高致使堤顶漫流破坏。这两种水流运动均为明流，与溃坝水流运动十分相似。因此，溃堤水流运动在数学上可以通过浅水方程(对一维问题即著名的圣维南方程)描述。堤身破坏形成溃口后，口门附近的河床将不断刷深形成冲刷坑，冲刷坑的发展将进一步影响堤脚的安全，并可能导致更大范围的堤防溃决。因此，准确地预测溃堤水流的行进过程及溃口下游冲刷坑的演化过程对抗洪抢险措施的决策实施具有指导意义。 随着计算机性能的提高以及数值计算技术的发展，对实际的溃决水流问题采用数值模拟与预报，不仅是可行的，而且是高效的。然而由于浅水方程的对流占优及退化的双曲型性质，尤其是溃口水流结构复杂，同时兼有急流与缓流，流场存在迅变或间断流动区域，如何较好的处理上述困难是模拟溃口流动的关键。迄今为止，关于堤身破坏形成溃口后口门附近冲刷坑发展过程的计算成果还很少见到。本文采用高分辨率格式的有限元法捕捉溃坝波的传播[3]，假定土体为无粘性沙，土体的冲刷输移满足泥沙连续方程和河床变形方程，建立了堤坝溃瞬间溃决后水流运动和冲坑发展的水动力学模型。给出了不同溃口宽度、不同流量和不同粒径土体组成条件下瞬时溃决的水流运动、冲坑发展及初步结论，为钢木土石组合坝的结构受力及稳定性分析提供水动力学条件。 2 基本方程 2.1 守恒形式的二维浅水方程 (1) (2) (3) 式中 qx、qy及u、v分别表示x、y方向的单宽流量和流速；λ=g/(C2h)；C为谢才阻力系数，由曼宁公式计算；水流涡粘性系数由vt=κu*h/6近似，其中κ为卡门常数，u*为摩阻流速；；水位函数ξ(x,y,t)由水深h(x,y,t)和床面底高程Zb(x,y,t)确定，即ξ(x,y,t)=h(x,y,t)+Zb(x,y,t)。 (4) (5) 式中 sk、s*k及ωk分别为第k粒径组悬移质泥沙的含沙量、挟沙力和沉速；φ=hsk；α为恢复饱和系数；泥沙紊动扩散系数εs假定与水流涡粘性系数相等；γ′s为淤积物干容重；Ns为悬沙的分组数。床沙起动临界流速采用张瑞瑾公式 Uck=(h/dk)0.14(17.6ρs-ρ/ρdk+6.05×10-710+h/dk0.72)1/2 (6) 式中 dk为第k粒径组泥沙的粒径；ρ、ρs为水、沙的密度。 2.3 水流挟沙力级配 水流中的泥沙来源于上游水流挟带和床沙紊动扩散进入，因此采用水流条件和床沙组成推求分组挟沙力的做法是较为合理的。本文采用李义天方法[4]计算挟沙力级配。 2.4 床沙级配计算 若已知各粒径组泥沙的冲淤厚度ΔHsk及总冲淤厚度ΔHs，则床沙级配的计算可分为以下两种情况： (1)ΔHs0，即发生淤积的情况，床沙活动层的级配由下式计算 ΔPbkt+Δt=ΔHsk+ΔPtbk(Htm-ΔHs)/Hmt+Δt (7) 式中 ΔPtbk、ΔPbkt+Δt分别为t、t+Δt时刻的床沙活动层级配，Htm、Hmt+Δt为相应时刻的床沙活动层的厚度。 (2)ΔHs0，即发生冲刷的情况，床沙活动层的级配由下式计算 ΔPbkt+Δt=ΔHsk+ΔPtbkHtm+｜ΔHs｜Δpremk/Hmt+Δt (8) (9) Mij =Cijqxj+Dijqxj+Fx (10) Mij=Cijqyj+Dijqyj+Fy (11) Mij=Cij(hs)j+Dij(hs)j+Fz (12) 式中 Mij为集中质量矩阵，Cij为对流矩阵，Dij为扩散矩阵，F为源