2019年6月
第46卷第11期 施工技术Construction Technology建筑技术开发Building Technology Development溃坝淹水及泥石流的数值模拟王福良，曹文明(水发规划设计有限公司，济南250014 )［摘 要］以有限差分法求解一含有源项的二维浅水波方程组，模拟中为了减少数值震荡，将数值通量以一阶精度与二阶精 度餉型态相结合，以符合二阶精度混和型态总变量消减(TVD)格式，搭配若干数值通量式，并加入干、湿交界面处理， 使其能有效处理在实际复杂地形底床坡度剧忍变化及干湿锋面的演变过程。［关键词］泥石流；浅水波方程组；总变量消减［中图分类号］TU201.4 ； U443.22 ［文献标志码］A ［文章编号］1001-523X (2019) 11-0070-02Numerical Simulation of Dam Breakwater and Debris FlowWang Fu-liang, Cao Wen-ming［Abstract ］ A two-dimensional shallow water wave equation system with source terms is solved by finite difference method. In order to reduce numerical oscillation, numerical flux is combined with first-order accuracy and second-order accuracy in order to conform to the second-order accuracy mixing state total variable subtraction (TVD) scheme. Several numerical flux restriction formulas are combined with the dry-wet interface treatment, which can effectively deal with the practical complexity. The sharp change of slope and the evolution of wet and dry front of shaped bed.［Keywords ］ low carbon concept ； interior art design ； application1研究目的与方法本研究主要以数值模拟的方法模拟流域中发生泥石流或 淹水时，其流动速度、堆积深度以及影响范围，利用二阶流变 模式配合有限差分法结合总变量消减(TVD)方法，建立二 维数值模式。利用有限差分法求解含有屈服强度的浅水波方程 式，引入若干单参数或双参数型函数如MUSCL, Superbee 等(表1),进行各种流场计算。最终研究成果希望能够为未 来泥石流等防治工程以及预警系统提供重要的依据。表1式分类内给定4个测点，用来测定当洪峰流量来时，其水深变化情形, 作为河川下游淹水潜势分析之用。测点坐标分别为G1 (0,240)、 G2 (480,1960)、G3 (2560,1960)与 G4 (39202560)。函数种类函数形式卩d=max {0,min［2,2rV+,0.5(r_+r+)］}MUSCLSuperbee^(r,r+)=max［0,mm(l ,2C,min(2F)］o4000 3 000 2000 1 000 02数值计算结果与讨论2.1案例讨论12.1.1计算条件概述在淹水模拟部分，模拟条件假定A水库受一突发大洪水 侵袭，水库溢淹，造成下游淹水。假设下游粗坑坝被冲毁，故 在模拟起始时将粗坑坝坝体移除，而模拟开始时粗坑坝内无 蓄水，渠道为干床状态。模拟范围为3960mx3960m的计算域, 包含不规则的山谷地形及蜿蜒的渠道。模拟初始条件曼宁系 数\"W0.3,模拟时间由0-61400S.由于入流条件采用变量流 模式给定，变量流入流量参考Hager (1985)所提出的变量流 方程式，水深％可由流量。算出：。©((护⑷)式中：Q为单宽流量；eP为单宽洪峰流量；tp为洪峰到 达时间。2.1.2结果讨论数值地形物理模型如图1所示，入流口位于地形右侧(A 水库坝体处)，出口位于左侧。地形由DIM资料构建，在地形收稿日期：2019-03-20作者简介：王福良(1986-),男，山东潍坊人，工程师，主要研究方向
为水工结构。
图1 A水库实际地形淹水流况的测点示意当时间为300s时，洪水波流入粗坑坝的山谷；700s时 洪水波流入粗坑坝蓄水区，其波峰有向四周扩散的流况；在 900s时其洪水已传达到蓄水区四周；之后，随着洪水传播, 流过粗坑坝体位置，假设无坝体存在，遂无溢淹坝体情形, 继续随着渠道向下游继续传播。图2为4个测点Gl, G2, G3, G4的水深变化历程图，历时 61400s。可以看出其最大洪峰水深各约为36m,13m,3.5m与 2m左右;随着时间的增加，约在18000s后水位逐渐趋于稳定, 持续保持恒定水位。而G2测点属于粗坑坝的蓄水区，所以水 位变化的恒定水位略高于入流的水位。2.2案例讨论22.2.1计算条件概述此案例中以B水库地区为主要仿真区域，进行溃坝淹水 模拟。DTM资料精确度为20mx20m,网格大小为102x102, 曼宁系数”值为0.095,给定的总流量为与前人比较分析，因 此参考利用FLO-2D模拟流量8000,总模拟时间为1800s,并 專立观测点对整体流况做进一步观察、分析，3个观测点的坐 标如下。P1 : (32, 67)oP2 : (36, 41)。• 70 •建筑技术开发Building Technology Development施工技术Construction Technology第46卷第11期2019年6月
时间图2模拟实际地形淹水下游各测点G1,G2,G3,G4
水深历程图(历时61400s)P3 : (33, 29)o观测点由坝体附近河道直至下游地区，依序选取作为观 测流速及水深的标准。2.2.2结果讨论首先利用软件Tecplot,将仿真结果以三维方式表现，由模 拟图中可以清楚看出模拟时间Z=720s时，水流开始向低地运动, 而模拟区域上方靠近坝体附近的山坡地，由于地势较高阻挡水 流运动，使得水流往图形的下方流动，虽然整体流况已经向下 游发展，但图中可明显看出坝体附近一山丘，由于地势较高因 此尚未被水流淹没”模拟时间Z=1440s时,下游河道因溃坝淹水而注入大量水体, 靠近坝体下方区域因疏流不及，已有部分水体向模拟区域的下 方移动。模拟时间Z=1800s时，虽然水流在高低差的影响之下迅 速地汇入河道，但由于河道疏流的速度不够迅速，因此无法 实时的排除过多的水量，淹没的范围也因此扩大。模拟时间 Z=3000~6000s时，整个模拟时间的后半段，其水体扩散的范围 虽有些许扩张但已相当的缓慢几乎停滞。3个观测点的水深历时分布图显示，其中Pi大约在700s时, 流体正好达到此处，因此水深迅速上升，而片的观测点较靠 近坝体，其整体水深分布在溃坝后先急速上升，并趋于稳定 后逐渐下降；而观测点P2与P3则是建立在河道附近，且两点 距离较近，因此在图形上较为相似，且水深分布也因为附近 流体的汇入而持续的升高，模拟时间6000s结束时，整体流况 与Pi相近，均在溃坝后先急速上升而后再逐渐下降。此软件(Tecplot)所绘制的图形虽然能以三维方式展现模 拟结果，但是仅能显示地形高低与基本流况，无法有效的表现 出模拟结果对于仿真区域的影响，因此在本文中利用ArcGIS 等相关软件将资料进行处理，以达到二维动力波与DIM资料 的数值模式结合的目的。由模拟图中可以发现，将模拟结果经过处理后，利用GIS 图层展现，除了更清楚的表示出整个仿真区域的地形外，还 能藉由仿真区域的各种建筑物作为淹没范围的标准，在模拟 开始后的12mm, B水库附近的主要道路已被淹没，且流体 已跨越河道往南方流动；在模拟时间Z=1440s以及/=1800s 时，整个仿真区域的下半部分已有多处遭到淹没，河道下 方许多工厂、民宅均已被淹没；整体二维的流况为其流况 已无剧烈的变化，由此可知流体正稳定的往边界流出。与 前人研究结果相比较，利用FLO-2D模拟流量8000n?/s, Z=1440s以及Z=1800s时的流况，以目前的模拟结果来看在下 游部分，其FLO-2D的模拟速度快于本文发展的模式结果，推 求其主要因素在于曼宁系数n值的影响。因为FLO-2D的曼宁 系数\"值会视计算的不稳定性，自行调整曼宁系数”值，所 以其曼宁系数n值为非定值的关系。本模式结合DTM资料所展现的模拟结果，除了在图形上
易于判别外，利用图形中各种地形地物以及流体各种物理条 件进行判断与分析，可为未来相关防治工程以及预警系统提 供重要的依据。2.3案例讨论32.3.1计算条件概述模拟条件延续案例1的后续发展，假定A水库受到一突发 大洪水侵袭，当水库溢淹发生的事件过后，土砂大量吸收水 分后产生泥石流。其余条件均与案例一相同，仅将溃坝淹水 相关物理参数改为泥石流，及模拟总时改为8000so 2.3.2结果讨论泥石流，或称为土石流，是指大量的松散土体与水体的 混合体，在重力作用下，沿河道或自然坡面由高处往低处流 动的现象。泥石流中的土体种类繁多，其颗粒大小的分布范 围很宽，从粘土至卵石甚至巨石。泥石流体绝大部份为水、 泥浆、砾石，大多发生在山区，其运动特性介于流体与固体 之间。泥石流中土体的体积浓度介于挟砂水流和滑动土体之间。 泥石流按照物质组成可分为泥流型泥石流、砾石型泥石流、 一般型泥石流；而本文中案例3主要以泥流型泥石流(又称泥 流)为模拟条件，是指泥流中土体物质主要由粘土、粉土和 砂所组成，很少砾石及卵石颗粒，其流速约在2~20m/s左右, 而泥石流速度最高约可达水流速度的10倍。当时间为250s时，泥流由入流点开始进入山谷，并流向 粗坑坝上游的蓄水区；时间为500s及700s时，泥流将粗坑坝 蓄水区注满后继续向下游运动，由此可以看出泥流在蓄水区 中满溢的状况，并向四周扩散开，且往附近一小型支流延伸； 在900s时泥流已到达到粗坑坝，并将粗坑坝注满；而1800s 时泥流到达出流口，此时可看出整个仿真区域均被泥流所注 满；而模拟结束8000s时，泥流大部分均已流出，仅粗坑坝 以下尚有泥流流动状态。在清水与泥流的比较上，由案例1清水流r=900s,与本 案例泥流?=900s对比，可观察出，泥流流速约为清水流速的 2倍多，与前述泥石流能达到10倍水流流速不同，此原因估计 与模拟地形的坡度有关。本文所釆用的地形坡度较平缓，影响整体流动速度的发 展，虽然流速上仅为清水流速的2倍，但仍比清水流流速快， 亦在范围值之内，尚属合理。3结论(1) 本模式不仅能处理清水流的问题，在泥流问题上亦 能加以分析处理，尤其在复杂的实际地形上进行模拟，许多 参数需经过相当多的参考资料与测试而订定，在模拟难度上 虽高，但仍可在合理范围内完成模拟。(2) 本文利用GIS等相关软件进行仿真结果资料处理, 并套用于模拟地形，其结果与一般绘图软件比较上，除可显 示地形、流况外，更可将原始地物、建筑、河道等显示出来， 在案例比较分析上优于一般的绘图软件。(3) 本模式与商业软件FLO-2D仿真结果比较，整体流况 均相当符合，使用上虽不如商用软件方便，但本文偏向学术 发展，且本模式可依地形以及各种条件进行修改，而FLO-2D 虽方便于工程上应用，但却无法修改其程序代码。参考文献[1] 刘洋，齐清兰，张力霆.尾矿库溃坝泥石流的演进过程及防护措施研究 [JJ.金属矿山,2015 (12) : 139-143.[2] 李书娜.尾矿坝溃坝模型及数值模拟研究[D],长沙:中南大学，2010.[3] 唐得胜.基于FLO-2D模型的不同频率泥石流数值模拟研究[D],成都：
成都理工大学，2014.[4] 马秋娟，唐阳，宿辉.泥石流启动过程试验与数值模拟研究[J].科学技
术与工程，2015, 15 (25) : 7-10.• 71 •
