1.3　堤坝溃决水流的特点及数值模拟方法

随着计算机技术及数值计算技术的进步，堤坝溃决水流的数值模拟技术得以迅速发展。目前堤坝溃决水流数学模型已成为洪水风险分析的重要工具，长期以来一直受到国内外学者的关注。如美国的DAMBREAK模型、荷兰的DELFT3D模型，均可对堤坝溃决水流进行模拟，丹麦的DHI也对堤坝溃决洪水问题进行了研究。近年来，随着计算机技术和空气动力学的飞速发展，很多比较新的计算方法被应用到堤坝溃决水流模拟中来，并成为水流数值计算领域的热点问题。

1.3.1　堤坝溃决水流的特点

堤坝溃决是突发的洪水灾害，下泄洪水向下游演进，具有陡峭间断的锋面，运动强度较大，造成下游水位的陡涨，上游因水库、河道水体的突然下泄，水位迅速下降，逆行波向上游传播，口门区域的流动具有较大的水头落差，在次临界流和超临界流间转换，属于非恒定含间断的自由表面流动。溃坝水流的流态分布如图1.3所示。

堤坝溃决水流一般同时兼有急流和缓流、涌波和稀疏波等流态，其洪水形态与一般的河道洪水相比，有以下特点（谭维炎，1998；Aureli等，2000）。

（1）洪峰流量大，水位高，紊乱漫滩的二维效应（斜激波、不规则水面），除了纵向水流外，横向扩散也很明显，因此除较窄的峡谷区外，一般按照二维问题处理。

（2）自然条件下的地形极不规则，容易形成流态过渡，急流、缓流、临界流常常同时发生，另外洪水流动过程中还常伴有涌波和水跃。

（3）在洪水演进过程中，水面梯度较大，特别是在口门附近常常存在着间断。

（4）堤坝溃决洪水一般都是在干床上演进，这给数值计算造成很大的难度，计算时动边界处理不当极易产生非物理振荡，造成计算失稳。

由于堤坝溃决水流的上述特点，通用的求解浅水方程组的差分方法常遭失败，其原因在于所用的格式通常不能同时满足虚假振荡的抑制和足够高的精度，不是过分耗散就是数值振荡。例如，Preissmann四点格式，在不做修正的情况下不能模拟跨临界流（Meselhe等，1997），Preissmann四点格式被用于FLDWAV模型和CHARIMA模型，在FLDWAV模型中使用了LPI（local partial inertia）方法来使得Preissmann四点格式在计算跨临界流时保持稳定，但它在动量方程中当Froude数接近1时忽略了对流项，在流动的非恒定性较弱时，LPI方法可以较精确地求解峰值流动，但在非恒定性很强时可能会引起很大的误差（Fread等，1998）。二维河道水流模拟的经典计算格式ADI在应用到堤坝溃决水流模拟时会出现数值振荡现象，为了消除数值振荡需要在控制方程中添加人工黏性项（Liang等，2006）。因此，在构建堤坝溃决水流洪水演进的格式时应具有如下要求：①不产生虚假的物理振荡现象，能够保持计算的稳定性；②能够精确地捕捉到水面间断现象，具有较高分辨率；③当时空步长趋向于无穷小时，数值解收敛于物理解。

1.3.2　数值模拟的主要途径及方法

1.3.2.1　主要途径

堤坝溃决水流的特殊性主要是指其所对应的物理流场中存在间断，该性质使得数值研究溃坝水流具有特定的困难，多数算法常常失效。正确模拟堤坝溃决水流这类强间断或大梯度流动的现象，主要有两种途径（Toro，2000；谭维炎，1998）：激波拟合法（shock－fitting method）和激波捕捉法（shock－capturing method）。在20世纪80年代以前，由于计算机条件的限制，大多采用激波拟合法模拟间断水面；80年代后，计算机技术的迅速发展为激波捕捉法模拟间断提供了条件。

激波拟合法的基本思想就是把间断看成是内部非连续的边界面，利用间断两侧的水力要素，在间断点直接引用代数形式的间断关系，拟合一个能根据此条件衔接两侧流动的间断，而在间断外的光滑区可以使用适合连续流的计算方法。该类方法精度较高，能够准确模拟间断位置和间断传播速度，主要缺点是计算太过复杂，编制程序非常不便，要不断追踪运动间断，同时该方法还要求所求的流体运动的流场结构已知，这在大多数情况下是困难的，因为流场事先未知。

激波捕捉法的基本出发点是采用计算方法所固有的数值耗散效应自动捕捉间断。若使用与守恒律微分方程组相容的守恒型差分格式，则所得数值解在间断两侧自动满足间断条件，因而不论解中是否存在间断，可以不加区别地统一进行计算，不必进行激波拟合的特殊处理。该类方法简单通用，容易实现，缺点是对间断位置与传播速度等间断关系的计算是近似的，与拟合法的优缺点恰好相反。

值得注意的是，对非守恒型的浅水方程使用激波捕捉法可能会导致某些错误（Toro，2000），如有些方法会在间断面附近产生非物理的伪振荡，不能反映数值解的真实情况（Liang等，2006）。Hou和Lefloch（1994）在理论上对错误的原因进行了探讨。他们指出，在计算间断时如果使用非守恒格式，格式即使收敛，也不是收敛到原方程的物理解。因此，为了正确利用激波捕捉法求解间断，应该使用守恒变量、守恒方程的形式和守恒型数值解法。

目前，自适应捕捉间断的数值解法主要有以下两种。

一种是把间断看成梯度很大的连续性特殊情况，这一方法的实现，可以在求解的微分方程组中人为地引进一定黏性作用的扩散项以平滑间断解，称为人工黏性法；或者在对微分方程进行离散时，有些格式本身就带有类似黏性的项，如Lax－Wendroff格式和Lax－Friedrichs格式，也可以使间断有一个光滑的过渡，称为格式黏性法，不过，该方法一阶精度格式会把间断的过渡区拉得过宽。

另一种是采用求解Riemann解的思想来完成（Godunov，1959）。Riemann思想的引入给数值间断解的模拟注入了新的活力，它不仅适用于光滑的古典解，更可以适应各种具有大梯度、大变形解的情况，该方法已经成为大梯度水面计算的主流方法之一。

1.3.2.2　主要方法

浅水动力学计算方法，按离散基本原理可以分为特征线法（MOC）、有限差分法（FDM）、有限元法（FEM）、有限体积法（FVM）。前3种方法在很多缓流问题的计算中取得了很大的成功，但不完全适合求解堤坝溃决水流等强间断的流动。有限体积法相当于守恒方程的直接离散，对由一个或多个控制体组成的任意区域，以至整个计算区域，都严格满足物理守恒律，不存在守恒误差，并且能正确计算间断。下面分别介绍这4种方法在计算间断方面的特点。

1.特征线法（MOC）

20世纪50年代初，林秉南提出一维水流计算的特征线法，迄今为止仍在不断改进。特征线法是求解双曲型偏微分方程的最精确的数值解法，其基本思想在于对一阶拟线性双曲型偏微分方程利用二维空间的特征理论，可导出两族特征曲面和相应的特征关系式，对特征关系式进行离散求解可得到变量的数值解（Hunt，1983；张家驹，1966）。

特征线法物理概念明确，数学分析严谨，计算精度较高，对于堤坝溃决水流这样不连续的现象，它的特征关系仍然存在，因此仍可以用特征线法求解，不过在间断处该方法不能直接计算间断，在间断点需采用激波拟合法使两侧衔接起来（韦春霞等，2003；陈景秋等，2004）。在数值求解时，采用差分法离散特征方程会带来守恒误差，当水流状态沿程变化较大时，非齐次项计算较繁，可能带来较大误差。

2.有限差分法（FDM）

有限差分法自20世纪50年代首次应用于模拟河道水流以来，至今仍是水动力学计算中应用最为广泛的方法（金旦华等，1965；胡四一等，1991）。差分法相当于点近似，是以泰勒级数展开为工具，对水流运动微分方程中的导数项用差分式来逼近，从而在每一个计算时段可以得到一个差分方程组。如果差分方程组可独立求解，则称为显格式，如果需联立求解，则称为隐格式。

目前工程上常用的差分方法有求解一维Saint－Venant方程组的Preissmann格式和求解二维浅水方程的交替方向隐（ADI）格式，前已叙述，由于这些算法都是基于连续性的数值逼近的思想，因此模拟强间断水面时均遭到失败。有些差分格式，由于本身带有平滑间断的项，所以可以计算间断水面，如Lax－Friedrichs格式、Lax－Wendroff格式等，但该类格式由于所含耗散太强，容易使得数值解过分光滑，使得间断在很宽的范围内被展平。MacCormack格式出现后，由于其耗散低、捕获激波能力强，并且具有二阶精度，因而很快取代了Lax－Wendroff格式，但是其不足之处是在激波附近数值解中有振荡现象产生，因此需要引入人工耗散项（汪迎春，2001）。

为了正确处理间断问题，很多差分格式采用了Riemann问题的思想。该类方法比较典型的是Godunov方法（Godunov，1959）。Godunov较好地考虑了物理流动特性，利用Riemann问题的精确解来计算网格边上的物理量，然而，这一方法在间断附近的解不够好且只有一阶精度，而且计算效率较低，Godunov同时证明了高于一阶精度的差分格式会在间断附近产生非物理的伪振荡。为了克服这一矛盾，Harten（1983a）基于微分守恒律首先证明了计算格式具有TVD（total variation diminishing）性质的充分条件并且具体构造了修正通量的高分辨率TVD格式，该格式能够很好地抑制数值振荡，体现了良好的激波捕捉能力，其不足之处是在极值点上只有一阶精度。Harten等（1986）提出了本质无振荡格式（essentially non－oscillatory，简称ENO格式）；Liu等（1994）提出了Weighted ENO格式（简称WENO格式）。这些都是更高阶的高分辨率差分方法，都可以较好地模拟间断。

采用Riemann问题思想的高分辨率差分方法可以很好地模拟间断形状，基本上解决了高阶精度差分格式在间断附近的非物理振荡问题，是数值求解的优良方法之一。但是，由于该类方法一般用于矩形网格或正交曲线网格，因此对于复杂边界的适应方面存在着不足。

3.有限元法（FEM）

有限元法从20世纪70年代开始应用于计算水力学中，该方法相当于体近似，其原理是分单元对解逼近，使微分方程空间积分的加权残差极小化。有限元法在数学上适于求解椭圆型方程组的边值问题，不适用于求解以对流为主的输运问题，同时有限元法缺乏足够耗散，捕捉锐利波形比较困难，不适用于计算间断。

针对普通有限元方法不能计算间断的特点，发展了间断有限元方法（discontinuous FEM）。特别是20世纪90年代以来，出现了计算间断能力较好的Runge－Kutta间断Galerkin方法（Cockburn等，1990，1998；Schwanenberg等，2004）。它既采用有限元弱解变分形式，又采用单元上的插值逼近，同时允许在时间和空间离散时存在间断。间断有限元法单元之间的连接更加精细和复杂，可以捕捉锐利波形，有效抑制虚假振荡，得到稳定的计算结果（Schwanenberg等，2004；李宏，2004；张修忠等，2001）。但是对非恒定流的计算，每一个时间步要求解一个大型的线性方程组，耗时较多，因而在工程中的应用受到很大限制。

4.有限体积法（FVM）

有限体积法是20世纪80年代以来发展起来的一种新型微分方程离散方法，综合了有限差分法和有限元法的优点。FVM与FEM一样，将计算区域划分成若干规则或不规则形状的单元或控制体，而不像经典的FDM那样要求网格有结构序列，因而处理过程具有较强的灵活性，能够满足处理复杂边界问题的需要，同时在间断问题的数值模拟方面显示了独特的效果（谭维炎，1999）。

FVM不是直接对方程组进行数值离散，而是从积分形式的守恒型方程组出发，采用非结构网格进行离散，在控制体边界上形成间断解的Riemann问题。在求出每个控制体边界沿法向的输入（出）的流量和动量通量后，对每个控制体分别进行水量和动量的平衡计算，得到计算时段末各控制体平均水深和流速。因此，FVM与FDM和FEM相比，其物理意义更为清晰。如果跨边界通量的计算只使用时段初值，则为显式FVM；反之，如果涉及时段始末的值，则为隐式FVM。因为跨控制体间界面输运的通量，对相邻控制体来说大小相等，方向相反，故对整个计算域而言，沿所有内部边界的通量相互抵消。对于由一个或多个控制体组成的任意区域，以至整个计算域，都严格满足物理守恒律，不存在守恒误差。

FVM的误差主要来自对界面数值通量的估算，目前已经发展了多种形式的求解界面数值通量的方法。最简单的方法为将界面处的数值通量看做相邻控制体形心处通量的平均，相当于守恒型差分中心格式；通量向量分裂格式（FVS）一般按特征符号进行分裂，以便在选择差分方向同时考虑信息传播方向，因而具有经典逆风格式的特点；利用Riemann解的Godunov格式是目前求解大梯度流动的主流格式，目前应用较多的Godunov格式有Roe格式（1981）、Osher格式（1982）、HLL格式（Harten等，1983b）等；另外，FDM中捕捉激波的高分辨率方法以及已经设计好的许多格式也可以直接利用到FVM中，如MacComack格式、QUICK系列格式等。

综上所述，FVM在求解大梯度水面的运动中具有独特的优势，从某种意义上说，现今的FVM体现了FEM的几何灵活性、MOC的精度和FDM的效率。因此，本书在构建模型时选择了基于Godunov格式的有限体积法。