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

堤防溃决的流动分析及冲刷坑计算张修忠，王光谦（清华大学水沙科学教育部重点实验室）摘要：以平面二维浅水方程和不平衡输沙理论作为建模基础，采用高分辨率有限元格式捕捉溃坝涌波的传播，建立了堤坝瞬间溃决后水流运动和河床冲刷的二维数学模型。计算了不同溃口宽度、不同流量和不同粒径土体组成条件下的堤坝瞬时溃决的水流演进和冲坑发展，为钢木土石组合坝的结构受力及稳定性分析提供水动力学条件。关键词：堤防溃决；冲刷坑；高分辨率有限元格式；数值模拟基金项目：国家自然科学基金及水利部联合资助重大项目，(项目编号：59890200)。作者简介：张修忠(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)。2.2不平衡输沙方程和河床变形方程(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)式中ΔPtbk、ΔPbkt+Δt、Htm、Hmt+Δt同前，ΔPremk为若干个记忆层内的床沙平均级配。2.5床沙级配的分层记忆模式为了模拟床沙的粗化和细化现象，模型将床沙划分为床沙活动层及其下面的分层记忆层两部分。床沙活动层的厚度为Hm，相应的级配为ΔPbk。分层记忆层可根据实际情况分n层，各层的厚度及相应的级配分别为ΔHn、ΔPnk。计算中，当河床发生淤积时，记忆层数增加，增加层的级配为t时刻的床沙活动层级配ΔPtbk。当河床发生冲刷时，根据冲刷量的大小，记忆层数相应减少，且级配作相应的调整。3基本方程的有限元离散及求解由于有限元法采用非结构化网格，便于处理复杂的实际问题，故本文对基本方程进行有限元离散，控制方程(1)～(4)的弱解形式经分部积分后可得如下的有限元空间半离散方程(9)Mij=Cijqxj+Dijqxj+Fx(10)Mij=Cijqyj+Dijqyj+Fy(11)Mij=Cij(hs)j+Dij(hs)j+Fz(12)式中Mij为集中质量矩阵，Cij为对流矩阵，Dij为扩散矩阵，F为源项，由下列各式计算Mij=∫ΩNiNjdΩ(13)Fx=-∫Ω(gh+λqx)widΩFy=-∫Ω(gh+λqy)widΩFz=-∫Ωαw(s-s*)widΩ式中N、w分别为插数值=N，则构成Galerkin法，Galerkin法等价于中心差分格式，因缺乏足够的耗散，对于对流占优问题，常导致数值振荡。为此，本文采用高分辨率格式对对流项进行重构，即通过引入几乎相等的扩散与反扩散以保证格式的高精度，同时利用限制因子保证影响系数的非负性及解的保单调性。为使计算收敛或加快收敛，离散中对源项F进行负坡线性化。离散后的方程(9)～(12)为常微分方程，可采用多种显式或隐式方法求解。本文对时间导数项采用C-N格式，对于离散的代数方程采用QMR[5]方法迭代求解，该方法具有节省内存，收敛快的优点。对离散后方程组顺序求解，即由对流输运方程(9)计算水深h，由二维对流扩散输运方程(10)、(11)和(12)计算u、v和s，最后由河床变形方程得到冲刷深度。4计算格式的检验取矩形平底无阻尼河道，底高程为0.0m，上游水位10m，下游水位1m，初时时刻水体静止，在x=0m处发生瞬间全溃。图1和图2给出了120秒时的水位与流速的数值解与stokes解析解[6]的对比关系。可以看出，本文的数值解既没有出现高频振荡，又避免了过大的数值耗散，与解析解吻合良好。图1水位对比关系Fig.1Comparisonofwatersurfaceelevations图2流速对比关系Fig.2Comparisonofvelocities图3溃坝计算域示意图Fig.3Computationaldomainofthedam-break图4冲深的时间的发展过程Fig.4Depthevolutionofthescourhole图5冲坑形状Fig.5Contourofthescourhole图6冲坑形状Fig.6Contourofthescourhole图7局部流速场Fig.7Localvelocityvector图8局部流速场Fig.8Localvelocityvector5溃决水流运动及冲坑发展的计算及结果分析本文计算了坝体溃决和堤防溃决两种情况。对于溃坝泄流，主要观察在上游水位不变、所有流动通过口门的情况下，口门附近的水流运动和冲坑发展情况。对于溃堤泄流，则主要观察在河道分流、口门上游水位变化的情况下的流动和冲刷情况。5.1溃坝泄流时流场及冲刷坑的计算及结果分析计算域如图3所示。x坐标：-1000～1500m,y坐标：-500～500m，坝厚30m。平底河床，糙率为0.015，土体中径为0.3mm。初始条件：水体静止，上游水深7.5m，下游水深0.5m。边界条件：上游给定水位7.5m，溃口及其上游闭边界给定无滑移条件，溃口下游给定流速的导数条件。计算了口门宽50m、100m和200m的三种工况。限于篇幅，只能给出部分结果图，更多的图形可查阅报告[7]。图4为口门宽50m时溃口中轴线的冲深随时间的发展过程；图5和图6为该口门宽度条件下600秒和7200秒时的冲深等势线图；图7和图8给出了口门宽100m条件下600秒和7200秒时刻的局部流速矢量场。为更好的比较不同口门宽度对冲坑演化的影响，便于找出冲深与口门宽度之间的关系，将不同口门宽度情况下冲坑的最大冲深及其所在位置制成表1。表1的第二行元素为溃口宽度，计算中的土体颗粒中径为0.3mm。表1不同溃口宽度的冲刷结果Table1Erosionresultsofdifferentwidthsofthebreach最大冲坑深度(m)冲坑最深点的位置(m)时刻(s)50100200501002001000.0910.160.1496.0120.0165.06000.6331.151.3096.8119.0165.012001.1812.012.5378.098.0142.018001.6402.673.5060.080.0120.036002.754.15.6045.080.099.072004.3756.145.062.0由表1可以看出：在上游水位相同的条件下，随溃口宽度的增加，冲坑深度增加；冲坑最深点距离溃口的位置也有所增加，该位置随时间发展逐渐靠近口门，最终稳定到某一位置。由冲深等势线图可以发现，冲坑在平面上呈椭圆形，且随着时间的发展逾趋明显。由表1及冲深随时间的发展过程图中还可以得出如下结论：冲坑的冲刷率先增加，而后减少，并最终趋于零，至冲坑不再发展；冲坑的上游边坡陡于其下游边坡，而且随着冲深的增加，冲坑的上游坡度有增大的趋势。5.2溃堤泄流的流场及冲刷坑计算及结果分析如图9所示，计算域包括上端河道、下端堤外区域及联结二者的口门三部分。河道宽800m，底坡为万分之二。堤宽为50m，口门及堤外均为平底，糙率取0.02。边界条件：河道进口给定流量过程，河道出口给定水位流量关系，该水位流量关系由谢才公式计算，溃口下游给定流速的导数条件。初始条件由河道水流计算确定。在相同的流量和土壤粒径条件下，计算了口门宽100m、200m、300m和500m四种情况的水流运动和冲坑发展；在保持溃口宽度不变的条件下，计算了不同流量和不同泥沙组成条件下的水流运动和冲坑发展。同样，受篇幅限制，只能给出部分结果图，更多的结果可查阅报告[7]。图10、图11给出了口门宽300m时，不同时刻的流速矢量场；图12、图13为口门宽300m时溃口中轴线的冲深随时间的发展过程；图14给出了溃口宽300m时溃口中轴线的流速随时间的变化关系。同样，为便于比较不同工况下冲坑的演化及探求某些规律，将不同工况的第7200s时的最大冲深以及其所在位置的计算结果汇总于表2与表3。图9溃堤计算域示意图Fig.9Computationaldomainofthedike-break图10T=100s的流速场Fig.10VelocityvectoratT=100s图11T=7200s的流速场Fig.11VelocityvectoratT=7200s图12冲坑的时间发展过程Fig.12Depthevolutionofthescourhole图13冲坑的时间发展过程Fig.13Depthevolutionofthescourhole图14溃口中轴线流速随时间过程Fig.14Velocityevolutionatthecentrallineofthebreach表2不同溃口宽度最大冲深Table2Maximumerosiondepthsofdifferentwidthofthebreach时间(s)