1.4　堤坝溃决水流数值模拟研究进展

1.4.1　溃口出流模型

溃口出流过程作为溃坝洪水演进数学模型的输入边界条件，其精度的高低直接关系到整个计算过程的成败。因此，在进行溃坝洪水计算时非常重要的一步是准确预测溃口的出流过程。但是由于溃口发展的影响因素众多（表1.1），溃决机理复杂，目前准确模拟溃口的发展过程还十分困难。

近几十年来，溃口出流过程的预测已经有了较多的成果。主要的方法总体来说可以分为以下两大类（Wahl，1998）。

一类为采用溃坝的历史数据资料得出的回归方程来计算出流过程，这类回归过程一般是建立溃口出流过程与库容和坝高的相关关系的表达式。Walder等（1997）通过已有溃坝资料分析认为，溃口的形状和尺寸与库区一样在溃坝过程中起着重要的作用，并且溃口一般都是典型的梯形断面，溃口宽度是深度的2～4倍。Kirkpatrick（1977）、Macdonald等（1984）、Costa（1985）、Froehlich（1995a，1995b）等都先后提出过相应的经验公式，Wahl（2001）曾对该类公式做过系统的总结。Xu等（2009）收集了8个国家和地区的182座土石坝的资料，其中75座资料非常齐全，根据这些数据采用回归分析的方法给出了相应的经验公式，这些公式分别适用于各种不同的溃决条件，可以计算溃口洪峰、峰现时间、最终溃口宽度等重要信息。值得一提的是，在使用这类方法时，一定要明确所采用的经验回归公式是否适用于所研究的水库大坝。

另一类为采用带有溃口发展物理机制的理论或半理论方法来预测出流过程。Cristofano（1965）建立了第一个模拟土石坝逐渐溃决过程的数学模型，模型假定溃口为梯形，不考虑边坡的侵蚀冲刷，溃口的底坡为常数，并等于坝体材料的内摩擦角。Harris和Wagner（1967）在Cristofano模型的基础上做了发展，不同的是他们假定溃口形状为抛物线形。Lou（1981）和Nogueira（1984）分别提出了各自的土石坝漫顶逐渐溃决的数学模型，模型最为显著的特点是对溃口形状不做假定，根据剪应力分析以及冲蚀机理确定溃口形状的解析解。Singh（1996）对大坝溃口的模拟技术进行了总结和归纳，并提出了自己的经验公式和有关溃口出流过程计算的数学模型。

随着溃口模拟技术的发展，国外已有很多成熟的溃口模型。这类模型基本都属于第二大类：在这其中又可分为两小类：一类为参数型溃口模型，另一类为物理－参数型溃口模型。所谓参数型溃口模型是指通过设定溃口尺寸、形状和发展时间，利用堰流公式来计算溃坝洪水出流过程（Manville，2001）。这些模型通常采用相对简单的参数化经验方法来估算溃口发展过程。DAMBRK就是这类模型的代表，1984年，美国国家气象局研发了溃坝模型DAMBRK（Fread，1984a，1988a），该模型的特点为所需要的参数少，最终的溃口尺寸由模型设定的确定溃口形状的参数和最终的溃口底宽决定，在溃决历时内，溃口形状呈线性发展，该模型没有泥沙输移计算，出流过程通过宽顶堰流公式计算。美国水文工程中心的HEC－1模型的溃口计算模块也是属于参数型的（Singh和Snorrason，1984）。

物理－参数型溃口模型指的是溃口的发展考虑了溃口冲刷的物理机制，通过泥沙输移公式来计算溃口侵蚀过程。这类模型相对较多，Breach模型就属于该类型，是当前世界上应用较广的土石坝溃决模型，同样是Fread研制的（Fread，1984b，1988b），模型可以模拟由漫顶和管涌引起的溃坝，坝体可以均质，也可以包含两种不同材料，分别构成心墙和外部区域。模型采用Meyer－peter＆Muller泥沙输移关系来计算泥沙输移量，该模型对内摩擦角和黏聚力的参数比较敏感。BEED模型同样属于物理－参数型，是由Singh和Scarlatos共同开发，模型中溃口被划分为两部分：一部分为坝顶的水平溃口段，另一部分为下游坡面上倾斜的溃口槽。断面假定为梯形，冲刷公式采用Einstein－Brown公式计算（Singh和Scarlatos，1989）。另外，BRDAM模型（Brown和Rogers，1981）、MIKE－11模型（DHI，1991）、SIMBA（Temple等，2005；Hanson等，2005）、HR－BREACH（Mohamed，2002）、FIREBIRD BREACH（Wang和Kahawita，2002）等均属于此种类型。

虽然已有很多学者致力于溃口出流数学模型的研究工作，但是由于对溃口机理的认识还不成熟，目前的溃口预测模型还存在着很大的不确定性，现有的预测水平下，溃口洪峰流量的预测误差在±50%范围内（Morris，2000），因此如何提高溃口模型预测结果准确度依然是该领域的一个难点问题。

1.4.2　溃决水流数值模拟

早期的堤坝溃决水流研究中，由于计算机水平的限制，溃决水流的计算条件通常经过很大的简化，然后推求其解析解，得出的解析解对理解溃坝洪水波演进的机理分析有很大的意义，但是这些技术无法应用于自然条件下具有真实地形的溃坝水流的模拟。进入20世纪80年代后，随着空气动力学和计算机技术的飞速发展，堤坝溃决水流数值模拟技术也开始焕发了新的生机和活力。

在当前的堤坝溃决水流数值模拟技术中，以求解Riemann近似解的思想为基础的Godunov格式逐渐成为主流（Godunov，1959）。该格式不仅适用于光滑的古典解，同时可以适应各种具有大梯度、大变形解的情况，能够精确地自动捕捉间断，由于该格式这种独特的优势，使其在溃坝水流一维、二维数值计算中得到了比较广泛的应用（史宏达等，2006；Toro，1999）。目前该方法已经成为大梯度水面数值模拟的主流方法，比较常用的Godunov格式有Roe格式、Osher格式和HLL格式等（Toro，1999）。

Glaister（1988）将在空气动力学中应用广泛的Roe格式应用于一维浅水方程的求解。Alcrudo等（1992）将一维的Roe格式推广到可以计算任意的河道大段面形状。Alcrudo和Garcia－Navarro（1993b）又将Roe格式应用到二维算例。Zhao等（1994）采用Osher格式建立了基于非结构网格的有限体积数学模型，并成功应用于Kissimmee河的水流模拟，随后Zhao等（1996）又将FDS格式应用于溃坝水流的模拟。胡四一等（1995）采用Osher格式，在非结构网格单元中引入了逆风的概念，成功模拟了长江口二维水流运动。Nugic（1995）指出在采用非结构网格离散守恒型浅水方程时，如果底坡项处理不好，模型在计算复杂地形上的静止水体时会产生虚假流动现象，如果计算中静水一直保持静止状态，则称计算格式满足C特性（又称和谐性）。该问题近年来成为浅水方程非结构离散求解的一个热点问题，很多学者对此进行了研究。如Zhou等（2001）提出了水面梯度法，潘存鸿等（2003）提出了水位方程法，Valiani等（2006）提出了底坡项分解法，Brufau等（2002）、王志力等（2005a）提出了底坡项特征分解法，这些方法代表了不同的研究思路，都取得了较好的效果。同时，一些学者把有限差分法的TVD概念引入有限体积法，构造了很多高效的有限体积法TVD格式（Burguete等，2001；Sanders，2001；Soares等，2002；Delis，2002，2003）。采用Godunov格式模拟堤坝溃决水流面临的另一个难题是干湿处理问题，天然实际地形下干床上的堤坝溃决水流模拟难度非常大（Toro，1999）。近年来，该问题已经成为溃坝水流研究中的热点问题。Valian等（2002）以有限体积法求解二维浅水波方程，采用HLL格式，应用于法国南部的Malpasset溃坝模拟，得到了较好的计算结果；Brufau等（2004）在Roe格式下采用有限体积法和零误差的干湿处理方法对Malpasset大坝的溃决水流进行了模拟；Liang等（2006）采用TVD－McCormack方法对Toce河物理模型的溃坝水流进行了计算；Liao等（2007）采用TVD－Lax－Wendroff格式对我国台湾的Cukeng坝的溃坝水流进行了模拟。Wang等（2013）采用了一种较新颖的干湿处理技术很好地解决了真实地形上的干湿边界处理问题。

综上所述，随着Godunov格式在水动力计算领域的发展应用，其面临的一些关键技术难题逐步得到了解决，使得该格式应用于工程实践的条件已基本成熟。

1.4.3　城区溃决水流数值模拟

随着社会经济的发展，城区人口以及财产分布更加集中，如果堤坝溃决洪水进入城市，将会给城市带来致命的损失。因此，近年来城区溃决水流数值模拟逐渐成为数值计算领域中的一个新的热点问题。

常见的引发城市洪水的原因有堤防的溃决或洪水漫顶，特大暴雨的发生以及城区排水设施的故障。对城市溃决洪水进行模拟时通常会遇到诸多困难，如城区街道上存在着汽车等阻水障碍物，社区和建筑物内部的蓄水容量以及城市的复杂几何形状等都难以精确量化，这与天然条件下的堤坝溃决水流有很大的不同。

如果所研究的区域是一个建筑物高度密集的城区，认为街道是最主要的水流通道，可以将计算区域概化成二维路网（Mignot等，2006a），建筑物密集城区的三维分布如图1.4所示。在街道上溃堤水流表现出很强的一维特性，而在十字路口和街道转弯处水流现象则比较复杂，因此一维模型的应用受到了很大的限制。Shettar等（1996）的研究认为，二维浅水模型能够正确模拟该处的水位和流场分布情况；Khan等（2000）则采用数学模型成功模拟了洪水在十字路口的分流现象；Nania等（2004）通过试验对坡度较陡的十字路口的分流现象进行了研究。

目前，不少学者采用二维浅水方程对密集城区的溃堤洪水进行了研究，如Zoppou等（1999）对澳大利亚Canberra市的供水水库的溃决淹没情况进行了模拟；Inoue等（2000）对日本Osaka城的溃堤洪水进行了研究；Haider等（2003）对1988发生在法国Nimes城的溃堤大洪水进行了模拟计算，取得了不错的结果；Calenda等（2003）对发生在Roma城的1870年的溃堤洪水进行了模拟；Mignot等（2006b）和Soares等（2006）等则通过城区物理模型试验准确验证了二维浅水方程对城区溃堤洪水进行模拟的可行性。国内目前还未见到这方面的报道。

1.4.4　一维、二维模型的耦合计算

河道一维水流模型由于其较高的计算效率和灵活性得到了广泛的应用，但是溃堤后的水流一般具有明显的二维性质，因此一维水流模型不再适用。如果直接采用二维模型来对河道水流和溃堤水流进行模拟，则需要耗费大量的计算机时，达不到实时预报洪水的目的。为此，将一维和二维水流模型进行耦合来模拟溃堤水流则可以发挥模型各自的优势。

与溃坝问题相比，溃堤问题的最大区别在于河道的侧向水流运动，这也是导致溃堤问题更加复杂的主要原因。从水流运动的角度上讲，溃堤出流运动与侧堰水流运动非常相似（张修忠，2003）。因此，对于溃堤洪水一维、二维的连接条件，应用最广泛的是堰流公式法。刘德平（1998）提到两种估算溃口流量过程的方法，即水文学方法和水力学方法，其中水力学方法是利用堰流公式进行计算的；冯民权等（2002）提出运用堰流公式作为边界条件对蓄滞洪区溃堤洪水演进进行数值模拟。

近几年，一维、二维耦合模型有了较大的发展。Dhondia等（2002）采用Sobek－Rural软件包对天然河流的漫顶洪水进行了模拟；李云等（2005）将一维、二维洪水嵌套模型成功应用于淮河临淮岗段；Lin等（2006）联合运用一维ISIS模型和二维DIVAST模型对伦敦附近的Greenwith港口的漫堤洪水进行了风险分析；Liang等（2007）采用动态链接技术将ISIS模型和DIVAST模型成功进行了链接，并通过算例验证了模型的质量守恒性，但是一维模型和二维模型的计算时间步长不能独立选取；Dushmanta等（2006）采用有限差分法的耦合模型对湄公河的漫顶洪水进行了成功模拟。