CN102436550B - 复杂边界及实际地形上溃坝洪水的自适应模拟方法 - Google Patents

Abstract

一种复杂边界及实际地形上溃坝洪水的自适应模拟方法，该方法包括对于变形后的偏微分方程的GRP格式和基于变分原理的自适应网格移动。前者通过在边界法向解一个广义黎曼问题，使得对于变形后的PDE的数值通量中包含了河床底部的影响，中心离散的源项也包含了流场的影响，在非结构三角形网格上精确保证静态流。后者实现网格位置的调整，在新网格上实现物理量的守恒性插值，不仅提高了非结构三角形网格对物理量的逼近程度，也能保证静态流。本发明开发的方法在模拟浅水波时能有效提高数值解的分辨率、计算精度和计算效率，可为溃坝洪水的水动力及富营养化过程模拟提供了一种新的技术方案。

Description

复杂边界及实际地形上溃坝洪水的自适应模拟方法

技术领域

本发明涉及到环境及水利工程技术领域，尤其是涉及一种复杂边界及实际地形上溃坝洪水的自适应模拟方法。

背景技术

我国病险水库多且分布广，其溃决洪水已成为国家防灾减灾工作中的一个突出问题。尤其是2011年的一号文件《中共中央国务院关于加快水利改革发展的决定》，进一步明确了新形势下水利的战略地位，其中强调到2020年基本建成防洪抗旱减灾体系，因此，为了分析溃坝洪水风险，实现致灾过程预报与灾害评估，在水利工程的设计和管理中溃坝洪水数值计算成为国内外重要研究课题之一。然而，溃坝洪水的流态较为复杂，通常包含激波，亚临界流，超临界流等区域且存在大梯度解，属于间断流动问题，模拟内容包括滞洪区和下游洪水波的传播速度、波形、波高、水位、流量等特征要素随时间变化的过程，目前主流的控制方程组采用浅水波方程。考虑到底部地形变化的影响，含源项的浅水波方程是一个时间相关的非线性双曲守恒律方程，要求溃坝洪水数值模拟格式具备和谐性，即保持流速为零、水位为常数的恒定状态。为了得到更高精度且考虑底部地形变化的影响，逐渐涌现出一些非常成功的格式，包括有限差分方法，有限元方法，间断Galerkin方法，空气动力学方法。比如，Li和Chen(2006)在结构网格下利用空气动力学的广义黎曼问题(GRP)格式求解浅水波方程，结合自由面梯度法(SGM)方法重构初始数据和维数分裂法处理二维问题，该格式实现精确保证和谐性质。然而，实际复杂边界上溃坝存在大量非凸区域，而传统的结构网格存在网格生成困难，难以实现复杂边界的贴边效果，有必要采用非结构三角形网格对区域进行剖分。但是在非结构网格下维数分裂遇到了困难，GRP方法难以直接推广，更没有证明在非结构网格下GRP格式是否同样满足和谐性及具有数值模拟的高准确度。

另一方面，溃坝洪水数值模拟求解主要是基于固定网格设计的格式。然而，溃坝洪水的水动力过程一般属于射流形式，与远离该区域的流场、浓度场在时间和空间都存在较大的差异，而我们往往又仅关心局部极值或变化大的梯度的洪水演进区域，显然采用粗糙的固定网格无法提高局部区域的计算分辨率，而采用一致的密网格会则会显著增加计算成本。因此，提出了自适应网格并应用到浅水波方程的数值模拟，其中，自适应加密网格方法是通过在局部加密网格来提高分辨率；自适应移动网格方法通过网格自适应调整来提高数值解的分辨率，并同时改善数值计算对计算资源的要求。目前自适应移动网格主要包括基于变分原理的移动网格方法，移动有限元方法，基于求解抛物型偏微分方程的移动网格方法，基于调和映射的移动网格方法等等。比如，Tang和Tang(2003)提出了结构网格上双曲守恒律方程的一般性自适应移动网格方法，该方法的网格移动方程简单且易于实现，并应用到了浅水波方程的计算。作为该方法学延伸，Chen等(2008)提出了非结构网格下多分量流的自适应移动网格方法。为了提高复杂边界及实际地形上溃坝洪水数值模拟的空间分辨率、计算精度和计算效率，有必要开发在非结构网格下自适应和谐性模拟方法。

发明内容

本发明的主要目的是提出一种在非结构三角形网格下的自适应和谐性模拟方法，为复杂边界及实际地形上溃坝洪水的水动力过程模拟提供新的模拟工具。

本发明的技术方案为一种复杂边界及实际地形上溃坝洪水的自适应模拟方法，将待模拟区域按网格进行处理，每个网格是一个三角形单元，通过守恒型的物理量U反映每个网格上水量信息，物理量U的向量形式如下，

$U=\left[\begin{array}{c}h+b\\ \mathrm{hu}\\ \mathrm{hv}\end{array}\right]$

其中，h表示水深，b表示河床底部高程，u表示坐标轴x方向的流速，v表示坐标轴y方向的流速；

按网格处理的具体实现包括以下步骤，

步骤1，根据待模拟区域，输入物理区域的初始网格和逻辑区域的初始网格；在时刻t_n，物理区域中第i个三角形单元E_i所在网格记为物理区域中第i个三角形单元E_i上物理量U的积分平均值记为令n＝0，计算在初始时刻t_n＝t₀的积分平均值

步骤2，进行移动网格，调整物理区域上的网格位置，移动任一网格的具体实现包括以下步骤，

步骤2.1，令迭代步数v＝0，v＝0时网格等于当前时刻t_n的网格v＝0时积分平均值等于当前时刻t_n的积分平均值

步骤2.2，计算表征物理量变化梯度大小的控制函数ω_i，计算公式如下，

其中算子表示在参考坐标下的梯度，参考坐标为逻辑区域的坐标；h表示水深，b表示河床底部高程，α和β为正参数；

步骤2.3，根据步骤2.2所得控制函数ω_i，通过以逻辑区域作为中间变量，移动网格至新的位置，新的位置记为网格

步骤2.4，在网格上更新迭代到第v+1步时的积分平均值

步骤2.5，判断是否满足预设的迭代结束条件，如果是则迭代结束，更新在下一时刻t_n+1的网格位置到网格更新当前时刻t_n的积分平均值转入步骤3；否则令迭代步数v＝v+1，转步骤2.3；

步骤3，根据步骤2的移动网格结果，求取下一时刻t_n+1的网格和积分平均值包括以下步骤，

步骤3.1，在网格上，以步骤2更新的积分平均值作为初始值，求取下一时刻t_n+1的积分平均值的近似解；

步骤3.2，如果下一时刻t_n+1＜TSTOP，TSTOP为预设的停止时刻，则令n＝n+1，以下一时刻为当前时刻返回执行步骤2，否则输出下一时刻t_n+1的网格所在位置和积分平均值作为模拟结果。

而且，步骤2.5中，采用或迭代步数v＝5作为预设的迭代结束条件。

而且，步骤2.3中，通过解离散的网格方程移动网格至新的位置，离散的网格方程如下，

其中，表示迭代的网格坐标位置的中间变量，表示迭代时的权重，μ_i为松弛因子，表示当网格x_i迭代到第v步时的邻居节点，N(i)表示网格x_i周围邻居节点的个数，且

其中，|ξ_ijξ_i|，|l_ij|分别表示逻辑区域上的网格边长及对偶网格边长，Δτ和σ是控制网格质量的变量，ω_ij是边l_ij的两个相邻单元上的控制函数w_i与w_j的平均值，|V_i|表示逻辑区域上以节点ξ_i为中心的对偶单元V_i的面积。

而且，步骤2.4中，通过守恒型插值公式得到积分平均值守恒型插值公式如下，

$\left|E_i^{[v+1]}\right|U_i^{[v+1]}=\left|E_i^{\left[v\right]}\right|U_i^{\left[v\right]}+\underset{j=1}{\overset{3}{\mathrm{\Σ}}}\hat{F}(U_{\mathrm{ij}}^L,U_{\mathrm{ij}}^R)$

其中，表示第v+1步时三角形单元E_i的测度，表示第v步时三角形单元E_i的测度，通量采用如下的迎风格式

$\hat{F}(U_{\mathrm{ij}}^L,U_{\mathrm{ij}}^R)=\max \left\{\right|D_{\mathrm{ij}}|,0\}{U}_{\mathrm{ij}}^L+\min \left\{\right|D_{\mathrm{ij}}|,0\}{U}_{\mathrm{ij}}^R$

其中，和是由单元平均值利用SGM方法重构出来的边l_ij的左右极限值，|D_ij|是网格移动到网格时边l_ij所扫过区域D_ij的测度，边l_ij表示三角形单元E_i的边。

而且，步骤3.1中，求取积分平均值的实现方式如下，

首先，以步骤2所得积分平均值作为初始值，利用自由面梯度法SGM重构分片线性的初始函数，初始函数见下式

其中，W(x，y)表示向量W的线性逼近函数，是向量W于当前时刻t_n在三角形单元E_i上的近似值，向量W＝(h+b，u，v)^T；r＝(x-x₀，y-y₀)表示位置坐标(x，y)相对三角形单元E_i中心(x₀，y₀)的相对坐标，梯度是由三角形单元E_i上的三个相邻单元E_i1、E_i2、E_i3上向量W的值，是一个选定的限制器；

然后，将在中点时刻边l_ij的中点处物理量代入下式，得到物理量U的源项S(U)的离散， $S\left(U\right)=\left[\begin{array}{c}0\\ {\mathrm{gb}(h+b)}_x\\ {\mathrm{gb}(h+b)}_y\end{array}\right],$ ((h+b)_x，(h+b)_y)表示自由面高程的梯度，

$S{\left(U\right)}_i^{n+\frac{1}{2}}={[0,{\mathrm{gb}}_i{(h_x+b_x)}_i^{n+\frac{1}{2}},{\mathrm{gb}}_i{(h_y+b_y)}_i^{n+\frac{1}{2}}]}^T$

其中，表示源项S(U)的近似值，g表示重力加速度，b_i表示在三角形单元E_i中心处河床底部高程；上标表示中点时刻Δt＝t_n+1-t_n表示时间步长；(b_x，b_y)表示河床底部高程b的梯度；h_x表示水深h的x方向变化率，h_y表示水深h在y方向的变化率；

最后，利用下式得到下一时刻t_n+1的积分平均值的近似解，

$U_i^{n+1}=U_i^n-\frac{\mathrm{\Δt}}{\left|E_i\right|}\underset{j=1}{\overset{3}{\mathrm{\Σ}}}{\mathrm{Fn}}_{\mathrm{ij}}{\left(U\right)}_{\mathrm{ij}}^{n+\frac{1}{2}}{\mathrm{\Δl}}_{\mathrm{ij}}+\mathrm{\ΔtS}{\left(U\right)}_i^{n+\frac{1}{2}}$

其中，Δt＝t_n+1-t_n表示时间步长，|E_i|表示三角形单元E_i的面积，Δl_ij为边l_ij的长度；和分别表示通量F(U)和与源项S(U)的近似值，通量 $F\left(U\right)=\left[\begin{array}{c}\mathrm{hu}\\ {\mathrm{hu}}^2+\frac{1}{2}{(h+b)}^2\\ \mathrm{huv}\end{array}\right].$

而且，在中点时刻边l_ij的中点处的物理量通过解广义黎曼问题得到近似值，具体方式如下，

$U_{\mathrm{ij}}^{n+\frac{1}{2}}=U_{\mathrm{ij}}^n+\frac{\mathrm{\Δt}}{2}{\left(\frac{\∂U}{\∂t}\right)}_{\mathrm{ij}}^n$

在当前时刻t_n边l_ij的中点处的物理量通过下式求出近似值

$U_{\mathrm{ij}}^n=\underset{t\→0^{+}}{\lim }U(0,t)=\mathrm{RM}\_\mathrm{Solver}(0;U_L,U_R)$

其中，RM_Solver(0；U_L，U_R)表示以U_L和U_R作为左右状态的黎曼问题解在时间轴上的值，U_L和U_R表示由SGM方法重构出来在单元边界的左右极限值；

时间方向导数通过以下关于水深时间变化率和水速时间变化率的二元一次方程组求出，

$\left\{\begin{array}{c}{a}_{-}^1{\left(\frac{\∂h}{\∂t}\right)}_{\mathrm{ij}}^n+a_{-}^2{\left(\frac{\∂u}{\∂t}\right)}_{\mathrm{ij}}^n=d_{-},\\ a_{+}^1{\left(\frac{\∂h}{\∂t}\right)}_{\mathrm{ij}}^n+a_{+}^2{\left(\frac{\∂u}{\∂t}\right)}_{\mathrm{ij}}^n=d_{+},\end{array}\right.$

其中，和d₊是常数。

本发明提供在非结构三角形网格下的自适应和谐性模拟方法，可以用于表征具有复杂边界及实际地形上溃坝洪水的水动力机理过程。并且具有如下优点：1)非结构三角形网格上的针对复杂边界及实际地形上溃坝洪水的广义黎曼问题(GRP)有限体积法具有高精度，而且能够精确的捕捉静态流，能够平衡复杂地形与流场的相互影响。2)针对溃坝洪水在非结构三角形网格自适应移动方法能够根据流场的特征，能够动态调整网格的空间分布，实现对激波，溃坝波和涌潮等局部极值或梯度变化大的重要区域的高分辨率捕捉。3)能够大大节约计算资源，也能够节约计算时间。本发明所提供非结构三角形网格上的高精度广义黎曼问题(GRP)有限体积离散格式具有特点：通过自由面梯度法进行初始数据重构，在单元边界中点处解一个广义黎曼问题，得到数值通量和源项的离散。由于GRP解中包含了河床地形的影响，从而实现了流场和河床地形的相互作用。又由于是对变形的浅水波方程进行离散，从而实现了对静态流的捕捉。本发明还提供了非结构三角形网格上的高精度广义黎曼问题(GRP)有限体积离散格式，通过自由面梯度法进行初始数据重构，在单元边界中点处解一个广义黎曼问题，得到数值通量和源项的离散。由于GRP解中包含了河床地形的影响，从而实现了流场和河床地形的相互作用；由于是对变形的浅水波方程进行离散，从而实现了对静态流的捕捉。本发明所提供的非结构三角形网格自适应移动方法通过解一个Euler-Lagrange方程来实现网格位置的调整，在新网格上基于几何方法实现物理量的守恒性插值，由于是对自由面高程在新网格上插值，从而该方法能够精确保证静态流。由于该方法将网格集中到人们关心的区域，从而节约计算资源，提高模拟效率的同时，提高了非结构三角形网格对物理量的逼近程度。

附图说明

图1为本发明中的非结构三角形网格示意图。

图2为本发明中的网格移动过程示意图。

图3为本发明应用实例二中平地二维溃坝在t＝7.2时网格以及自由面高程h+b的等值线示意图，其中图3a为固定的粗网格，图3b为固定的密网格，图3c为移动的粗网格，图3d为在粗网格下的计算结果，图3e为在固定的密网格下的计算结果，图3f为在移动的粗网格下的计算结果。

图4为本发明应用实例二中平地二维溃坝在t＝7.2时在横断面y＝132m处的水流速度示意图；

图5为本发明应用实例三中圆形溃坝在是t＝2时网格以及自由面高程h+b的等值线示意图，其中图5a为固定的粗网格，图5b为固定的密网格，图5c为移动的粗网格，图5d为在粗网格下的计算结果，图5e为在固定的密网格下的计算结果，图5f为在移动的粗网格下的计算结果。

图6为本发明应用实例三中圆形二维溃坝在t＝2.0时在横断面y＝25m处的自由水面高程示意图。

图7为本发明应用实例三中圆形维溃坝在t＝0，0.4，0.8，1.2，1.6时的网格变化示意图，其中图7a为t＝0时的网格图，图7b是t＝0时相应的自由面高程h+b的等值线示意图，图7c为t＝0.4时的网格图，图7d是t＝0.4时相应的自由面高程h+b的等值线示意图，图7e为t＝0.8时的网格图，图7f是t＝0.8时相应的自由面高程h+b的等值线示意图，图7g为t＝1.2时的网格图，图7h是t＝1.2时相应的自由面高程h+b的等值线示意图，图7i为t＝1.6时的网格图，图7j是t＝1.6时相应的自由面高程h+b的等值线示意图。

具体实施方式

为便于理解本发明技术方案的实施和效果起见，以下结合附图提供实施例的具体方案说明：

将待模拟区域按网格进行处理，每个网格是一个三角形单元，通过守恒型的物理量U反映每个网格上水量信息，物理量U的向量形式如下，

$U=\left[\begin{array}{c}h+b\\ \mathrm{hu}\\ \mathrm{hv}\end{array}\right]$

其中，h表示水深，b表示河床底部高程，u表示坐标轴x方向的流速，v表示坐标轴y方向的流速。

按网格处理的具体实现包括以下步骤，具体实施时可采用计算机软件技术实现自动处理：

步骤1，输入待模拟区域的物理区域的初始网格和逻辑区域的初始网格；在时刻t_n，物理区域中第i个三角形单元所在网格记为物理区域中第i个三角形单元上物理量U的积分平均值记为令n＝0，计算在初始时刻t_n＝t₀的积分平均值

物理区域的初始网格划分可以通过现有软件实现，把整个待模拟区域分成多个三角形单元。逻辑区域是为了数值计算构造出来的一个虚拟区域，目的是为控制物理区域上网格的空间分布。他本身只具有数学意义，而没有物理意义。逻辑区域的初始网格取和物理区域的初始网格一样即可。

步骤2，进行移动网格，调整物理区域上的网格位置。本发明设定每次要更新物理量的时候都调整一下网格的位置。实施例中，移动任一网格的具体实现包括以下步骤：

步骤2.1，令迭代步数v＝0，v＝0时网格等于当前时刻t_n的网格v＝0时积分平均值等于当前时刻t_n的积分平均值

步骤2.2，计算表征物理量变化梯度大小的控制函数ω_i，计算公式如下，

其中算子表示在参考坐标下的梯度，h表示水深，b表示河床高程，α和β为依赖于具体问题的正参数。这里的参考坐标是逻辑区域的坐标系统。

步骤2.3，根据步骤2.2所得控制函数ω_i，通过以逻辑区域作为中间变量，移动网格至新的位置，新的位置记为网格实施例中通过解离散的网格方程移动网格至新的位置，参见式(17)。

步骤2.4，在网格上更新迭代到第v+1步时的积分平均值实施例中通过守恒型插值公式得到积分平均值守恒型插值公式见式(22)。

步骤2.5，判断是否满足预设的迭代结束条件，如果是则迭代结束，更新在下一时刻t_n+1的网格位置到位置更新当前时刻t_n的积分平均值转入步骤3；否则令迭代步数v＝v+1，转步骤2.3。

步骤3，根据步骤2的移动网格结果，求取下一时刻t_n+1的网格和积分平均值包括以下步骤：

步骤3.1，在网格上，以步骤2更新的积分平均值作为初始值，求取下一时刻t_n+1的积分平均值的近似解。实施例求取积分平均值的实现方式为，以步骤2所得积分平均值作为初始值，利用自由面梯度法SGM重构分片线性的初始函数，参见式(9)；然后将广义黎曼问题的解代入式(7)得到源项的离散，广义黎曼问题的解见式(14)和式(15)；然后再利用式(5)更新物理量。

步骤3.2，如果下一时刻t_n+1＜TSTOP，TSTOP为预设的停止时刻，则令n＝n+1，以下一时刻为当前时刻返回执行步骤2，否则输出下一时刻t_n+1的网格位置和积分平均值作为最后结果。

为便于实施参考起见，本实施例的具体计算实现详细说明提供如下：

1.溃坝洪水问题的数学模型描述

本发明的主要目的是提出一种在非结构三角形网格下的自适应数值模拟方法，为复杂边界及实际地形上溃坝洪水的水动力过程模拟提供新的计算工具。并重点分析流场和复杂地形之间的和谐效应。因此仅考虑只有底部地形变化的影响的二维浅水波方程，该方程可以用如下偏微分方程描述：

U_t+F(U)_x+G(U)_y＝S(U)    (1)

其中，

$U=\left[\begin{array}{c}h\\ \mathrm{hu}\\ \mathrm{hv}\end{array}\right],$ $F\left(U\right)=\left[\begin{array}{c}\mathrm{hu}\\ {\mathrm{hu}}^2+\frac{1}{2}g{h}^2\\ \mathrm{huv}\end{array}\right],$ $G\left(U\right)=\left[\begin{array}{c}\\ \mathrm{hv}\\ \mathrm{huv}\\ {\mathrm{hv}}^2+\frac{1}{2}{\mathrm{gh}}^2\end{array}\right],$ $S\left(U\right)=\left[\begin{array}{c}0\\ -{\mathrm{ghb}}_x\\ {-\mathrm{ghb}}_y\end{array}\right],---\left(2\right)$

式中，t表示时间，h表示水深，x、y分别表示两个坐标轴方向，u、v分别表示x、y方向的流速，g表示重力加速度，b表示河床底部高程，(b_x，b_y)表示河床底部高程的梯度。U为守恒型物理量，F(U)、G(U)分别表示在x、y方向的通量，S(U)为源项。

U_t表示物理量U的时间变化率，F(U)_x表示x方向通量在该方向的空间变化率，G(U)_y表示y方向通量在该方向的空间变化率。

此处忽略了水应力、底部摩阻力、科氏力和粘性力等的作用。为了便于数值格式的设计，本发明采用(1)的一种等价的微分方程：

$U=\left[\begin{array}{c}h+b\\ \mathrm{hu}\\ \mathrm{hv}\end{array}\right],$ $F\left(U\right)=\left[\begin{array}{c}\mathrm{hu}\\ {\mathrm{hu}}^2+\frac{1}{2}{(h+b)}^2\\ \mathrm{huv}\end{array}\right],$ $G\left(U\right)=\left[\begin{array}{c}\\ \mathrm{hv}\\ \mathrm{huv}\\ {\mathrm{hv}}^2+\frac{1}{2}{g(h+b)}^2\end{array}\right],$ $S\left(U\right)=\left[\begin{array}{c}0\\ {\mathrm{gb}(h+b)}_x\\ {\mathrm{gb}(h+b)}_y\end{array}\right],---\left(3\right)$

式中，h+b表示自由面高程，((h+b)_x，(h+b)_y)表示自由面高程的梯度。

2.非结构三角形网格上的广义黎曼问题(GRP)有限体积法

本发明的第一部分是设计对方程(1)、(3)在非结构三角形网格上的高精度有限体积离散，该格式是以控制体为中心的格式，边界处的数值通量采用GRP格式来计算，该数值通量考虑了地形的影响，GRP问题的解反过来用于底部地形的逼近，从而实现复杂河床地形与流场的相互作用和平衡。

设三角形E_i为二维空间R²里的多边形物理区域的第i个单元，R²表示时域。该单元的第j个顶点记作x_ij，与这个顶点相对的边记作l_ij，以l_ij作为公共边的相邻单元记为E_ij，n_ij表示l_ij的单位法向量，方向从E_i指向E_ij。如图1，三角形单元E_i的第1个顶点记作x_i1，与这个顶点x_i1相对的边记作l_i1，以l_i1作为公共边的相邻单元记为E_i1，n_i1表示l_i1的单位法向量，方向从E_i指向E_i1；三角形单元E_i的第2个顶点记作x_i2，与这个顶点x_i2相对的边记作l_i2，以l_i2作为公共边的相邻单元记为E_i2，n_i2表示l_i2的单位法向量，方向从E_i指向E_i2；三角形单元E_i的第3个顶点记作x_i3，与这个顶点x_i3相对的边记作l_i3，以l_i3作为公共边的相邻单元记为E_i3，n_i3表示l_i3的单位法向量，方向从E_i指向E_i3。

在时间t＝t_n时，物理量U在三角形单元E_i的积分平均值记作河床底部高程b用分片线性的全局连续函数来近似，设河床底部高程b在三角形单元E_i的顶点x_ij的值为b(x_ij)，则在三角形单元E_i内及单元边界上，河床底部高程b可用线性函数插值得到。

方程(1)在三角形单元E_i上积分可得，

$\left|E_i\right|\frac{d}{\mathrm{dt}}{U}_i=-\underset{j=1}{\overset{3}{\mathrm{\Σ}}}{\∫}_{l_{\mathrm{ij}}}{\mathrm{Fn}}_{\mathrm{ij}}\left(U\right)\mathrm{ds}+{\∫}_{E_i}S\left(U\right)\mathrm{dx}---\left(4\right)$

其中，|E_i|表示三角形单元E_i的面积，ds表示单元边界的长度微元，表示单元边界的单位外法向量，表示在单位外法向量n_ij的数值通量。边上的积分和面上的积分都采用中点积分，时间方向导数采用中心差分。因此，在新时刻t_n的物理量U在三角形单元E_i的积分平均值更新为可采用如下形式

$U_i^{n+1}=U_i^n-\frac{\mathrm{\Δt}}{\left|E_i\right|}\underset{j=1}{\overset{3}{\mathrm{\Σ}}}{\mathrm{Fn}}_{\mathrm{ij}}{\left(U\right)}_{\mathrm{ij}}^{n+\frac{1}{2}}{\mathrm{\Δl}}_{\mathrm{ij}}+\mathrm{\ΔtS}{\left(U\right)}_i^{n+\frac{1}{2}}---\left(5\right)$

其中，Δt＝t_n+1-t_n表示时间步长，Δl_ij为边l_ij的长度；和分别表示通量F(U)和与源项S(U)的近似值，可以通过下面的公式求解，

${\mathrm{Fn}}_{\mathrm{ij}}{\left(U\right)}_{\mathrm{ij}}^{n+\frac{1}{2}}={\mathrm{Fn}}_{\mathrm{ij}}\left(U_{\mathrm{ij}}^{n+\frac{1}{2}}\right),$ 其中 $U_{\mathrm{ij}}^{n+\frac{1}{2}}=U(x_{l_{\mathrm{ij}}},t_n+\frac{\mathrm{\Δt}}{2})---\left(6\right)$

为边l_ij的中点，表示物理量U在中点处及中点时刻的解，即在中点时刻时，边l_ij的中点处的物理量。

$S{\left(U\right)}_i^{n+\frac{1}{2}}={[0,{\mathrm{gb}}_i{(h_x+b_x)}_i^{n+\frac{1}{2}},{\mathrm{gb}}_i{(h_y+b_y)}_i^{n+\frac{1}{2}}]}^T---\left(7\right)$

其中b_i，可用下面的中心格式离散，

$\left\{\begin{array}{c}{b}_i=\frac{1}{3}\underset{j=1}{\overset{3}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}\\ {(h_x+b_x)}_i^{n+\frac{1}{2}}=\frac{1}{\left|E_i\right|}\underset{j=1}{\overset{3}{\mathrm{\Σ}}}(h_{\mathrm{ij}}^{n+\frac{1}{2}}+b_{\mathrm{ij}})n_{\mathrm{ij}}^x{\mathrm{\Δl}}_{\mathrm{ij}}\\ {(h_y+b_y)}_i^{n+\frac{1}{2}}=\frac{1}{\left|E_i\right|}\underset{j=1}{\overset{3}{\mathrm{\Σ}}}(h_{\mathrm{ij}}^{n+\frac{1}{2}}+b_{\mathrm{ij}})n_{\mathrm{ij}}^y{\mathrm{\Δl}}_{\mathrm{ij}}\end{array}\right.---\left(8\right)$

其中，表示单元边界的单位法向量n_ij在x方向的分量，表示单元边界的单位法向量n_ij在y方向的分量；b_i表示在三角形单元E_i中心处河床底部高程，b_ij表示单元边界中点处河床底部高程。(b_x，b_y)表示河床底部高程b的梯度；h_x表示水深h的x方向变化率，h_y表示水深h在y方向的变化率。上标标识中点时刻Δt＝t_n+1-t_n表示时间步长。

因此，估算和的重点是求出在中点时刻的近似解本发明采用GRP格式来求解，其解法可以分成如下三个步骤来完成：

第一步采用自由面梯度法(SGM)对向量W＝(h+b，u，v)^T进行初始数据重构：

其中，上标符号‘T’表示的是向量的转置，W(x，y)表示向量W的线性逼近函数，是向量W于t_n时刻在三角形单元E_i上的近似值，r＝(x-x₀，y-y₀)表示位置坐标(x，y)相对三角形单元E_i中心(x₀，y₀)的相对坐标，梯度是由三角形单元E_i上的三个相邻单元E_i1、E_i2、E_i3上W值，利用现有的格林公式计算得到。是一个选定的限制器，当恒等于0时，物理量在控制体内离散为分片常数，这时候格式空间只具有一阶精度；当恒等于1时，格式具有二价精度，但是在激波等强间断附近，将出现数值扰动，需要采用斜率限制器来控制这种不稳定性。

实施例采用迎风型限制器，即：

其中，和表示向量W在三角形单元E_j的相邻单元的局部极大值和极小值，

$W_j^{\max }=\underset{k}{\max }(W_j^n,W_{\mathrm{jk}}^n),$ $W_j^{\min }=\underset{k}{\min }(W_j^n,W_{\mathrm{jk}}^n),---\left(11\right)$

其中，k取遍三角形单元E_j的相邻三角形单元。表示t_n时刻向量W在三角形单元E_j的相邻单元E_jk的近似值，表示t_n时刻向量W在三角形单元E_j的近似值。

第二步在三角形单元E_i的单元边界的外法向解一个局部的广义黎曼问题(GRP)：

外法向的一维方程有如下形式，

$\frac{\∂U}{\∂t}+\frac{\∂F\left(U\right)}{\∂\mathrm{\ξ}}=S\left(U\right)={(0,\mathrm{gb}{h}_x,0)}^T---\left(12\right)$

其中，上标符号‘T’表示的是向量的转置，x表示沿外法方向的局部坐标。

初始数据是利用式(9)重构出来的一个在外法线方向的线性函数

$W(\mathrm{\&xi;},0)=\left\{\begin{array}{cc}{W}_L+\mathrm{\&xi;}{\left(\frac{\&PartialD;W}{\&PartialD;\mathrm{\&xi;}}\right)}_L,& \mathrm{\&xi;}<0,\\ W_R+\mathrm{\&xi;}{\left(\frac{\&PartialD;W}{\&PartialD;\mathrm{\&xi;}}\right)}_R,& \mathrm{\&xi;}>0.\end{array}\right.---\left(13\right)$

其中，W_L和W_R分别表示W(x，y)在边界l_ij的中点处的左右极限值。和分别表示向量W在外法方向变化率在边界l_ij的中点处的左右极限值。

解式(12)和(13)可以得到具体通过解广义黎曼问题得到的实现方式如下：

其中，在当前时刻t_n边l_ij的中点处的物理量可以通过标准的齐次方程的黎曼问题解获得近似值，

$U_{\mathrm{ij}}^n=\underset{t\&RightArrow;0^{+}}{\lim }U(0,t)=\mathrm{RM}\_\mathrm{Solver}(0;U_L,U_R)---\left(14\right)$

其中RM_Solver(0；U_L，U_R)表示的是以U_L和U_R作为左右状态的黎曼问题解在时间轴(t轴)上的值，U_L和U_R表示由SGM方法重构出来在单元边界的左右极限值，SGM方法见式(9)。

时间方向导数可以通过利用解的结构以及初始条件解析的构造出来，只需构造一个关于水深时间变化率和水速时间变化率的二元一次方程组，

$\left\{\begin{array}{c}{a}_{-}^1{\left(\frac{\&PartialD;h}{\&PartialD;t}\right)}_{\mathrm{ij}}^n+a_{-}^2{\left(\frac{\&PartialD;u}{\&PartialD;t}\right)}_{\mathrm{ij}}^n=d_{-},\\ a_{+}^1{\left(\frac{\&PartialD;h}{\&PartialD;t}\right)}_{\mathrm{ij}}^n+a_{+}^2{\left(\frac{\&PartialD;u}{\&PartialD;t}\right)}_{\mathrm{ij}}^n=d_{+},\end{array}\right.---\left(15\right)$

其中，和d_-是只依赖W_L，b_L和的常数，b_L和表示沿边界外法向的左侧的河床底部高度及其空间变化率。和d₊是只依赖W_R，b_R和的常数，b_R和表示沿边界外法向的左侧的河床底部高度及其空间变化率。然后再计算具体的计算步骤参考可参考现有文献：Li JQ and G.X.Chen，The generalized Riemann problemmethod for the shallow water equations with bottom topography.International Journal forNumerical Methods in Engineering.65(2006)834-862.

第三步将计算出来的单元边界中点处在中点时刻的物理量的近似值代入式(7)得到源项S(U)的近似值的离散，然后再利用式(5)更新物理量。

以下分析本发明所设计的非结构三角形网上的GRP格式的能够精确保证静态流。不妨令时间t＝t_n时，三角形单元E_i所在河床底部x、y方向的水流速度自由面高程 由SGM方法知初始数据重构出来的多项式(9)满足又由河床高程函数的线性离散知，式(13)中的W_L＝W_R， 由广义黎曼问题方法知且代入式(5)可知时间t＝t_n+1时，三角形单元E_i所在河床底部x、y方向的自由面高程即在非结构三角形网上的GRP格式精确保证静态流的模拟。

3.非结构三角形网格上的自适应移动网格方法

本发明的第二个重点是根据溃坝洪水的流场特征，动态调整网格的空间分布，实现对激波，溃坝波和涌潮等局部极值或梯度变化大的重要区域的高分辨率捕捉。

本发明采用的是基于变分原理的网格移动方法，该方法通过解一个能够控制网格位置的网格偏微分方程来调整网格位置，再将物理量在新网格上插值。该方法能够保证精确捕捉静态流，并且能够提高数值解的分辨率，节约计算时间和计算资源。具体方法可参考：Tang HZ，TTang，Adaptive mesh methods for one-and two-dimensional hyperbolic conservation laws.SIAMJournal on Numerical Analysis.41(2003)487-515.及Chen GX，HZ Tang，PW Zhang.Second-orderaccurate Godunov scheme for multicomponent flows on moving triangular meshes.Journal of ScientificComputing.34(2008)64-86.

3.1基于变分原理的网格移动

逻辑区域取为直角坐标系ξ下区域Ω_i，并且逻辑区域Ω_i上的网格划分T_i具有与物理区域Ω_p的网格划分T相同的数据结构。从逻辑区域Ω_i到物理区域Ω_p的映射x＝x(ξ)，ξ∈Ω_i由变分原理得到网格移动方程，即Euler-Lagrange方程所决定，方程如下：

其中，表示逻辑区域上的梯度算子。在本发明中，表示控制函数，用来表征物理量变化大小的程度，进而控制网格的疏密程度。其中v表示的是网格迭代步数，表示第v步所得物理量U在三角形单元E_i的积分平均值。

首先考虑区域内部的网格移动，对于区域边界上的网格的同步移动将在后面说明，这里首先假设区域边界上的网格是固定不动的，对网格方程(16)的离散采用有限体积格式，离散后得到的非线性方程采用松弛Jacobi迭代法求解：

其中，表示Jacobi迭代的网格坐标位置的中间变量，表示迭代时的权重，μ_i为松弛因子。表示当网格迭代到第v步时节点x_i的邻居节点，N(i)表示该节点周围相邻节点的个数，且

其中，|ξ_ijξ_i|，|l_ij|分别表示逻辑区域上的网格边长及对偶网格边长，Δτ和σ用来控制网格的质量，在本发明取Δτ＝0.5和σ＝0.3。ω_ij是边l_ij的两个相邻单元上的控制函数w_i与w_j的平均值，控制函数根据步骤2.2中计算公式得到，|V_i|表示逻辑区域上以节点ξ_i为中心的对偶单元的面积。又由于网格的移动的目的是为了调整网格来提高网格对解的逼近程度，所以并不需要精确求解，于是这里以网格移动量，即第v+1步与第v步网格位置的差或者迭代步数v＝5作为迭代结束的标准。

由于间断或者人们关心的区域也会移动到物理区域边界上，所以物理区域边界上的网格点也需要和区域内部的网格同时移动。这里仅考虑的多边形物理区域，则边界线是一些直线段，于是把(16)投影到这些线段上面得到是一个一维的网格方程，如果这条线段平行与x轴，则(16)的第一个分量方程便是我们要离散的网格方程。

3.2基于几何的守恒性插值或重映

通过解网格移动方程(16)，将网格从旧的位置移动到了新的位置，这时需要将旧网格上的解的积分平均值插值或重映到新网格上，从而得到在新网格上的守恒量平均值而且这种插值满足如下的离散的全局守恒意义，

$\underset{i}{\mathrm{\&Sigma;}}\left|E_i^{\left[v\right]}\right|U_i^{\left[v\right]}=\underset{i}{\mathrm{\&Sigma;}}\left|E_i^{[v+1]}\right|U_i^{[v+1]}---\left(19\right)$

即在网格迭代第v步时的总质量与网格迭代第v+1步时的总质量相等，其中i遍历所有的三角形单元；上式中表示在第v步的三角形单元E_i，表示在第v+1步的三角形单元E_i。参见图2，第v步时，三角形单元E_i的第1个顶点在第2个顶点在第3个顶点在第v+1步时，三角形单元E_i的第1个顶点在第2个顶点在第3个顶点在

首先计算网格移动到网格时，边l_ij所扫过区域D_ij的测度|D_ij|，即带符号的面积。不失一般性，第一条边l_i1扫过的区域D_i1的测度|D_i1|可以通过下式计算

$\left|D_{i1}\right|=\frac{1}{2}((x_{i3}^{[v+1]}-x_{i2}^{\left[v\right]})(y_{i3}^{\left[v\right]}-y_{i2}^{[v+1]})-(y_{i3}^{[v+1]}-y_{i2}^{\left[v\right]})(x_{i3}^{\left[v\right]}-x_{i2}^{[v+1]}))---\left(20\right)$

其中，下标i2、i3分别表示边l_i1的两个顶点编号，表示第v步时顶点的两个坐标分量，表示第v步时顶点的两个坐标分量，表示第v+1步时顶点的两个坐标分量，表示第v+1步时顶点的两个坐标分量。第二条边l_i2扫过的区域D_i2的测度|D_i2|，和第三条边l_i3扫过的区域D_i3的测度|D_i3|同理计算即可。该测度计算具有如下性质

$\left|E_i^{[v+1]}\right|=\left|E_i^{\left[v\right]}\right|+\underset{j=1}{\overset{3}{\mathrm{\&Sigma;}}}\left|D_{\mathrm{ij}}\right|---\left(21\right)$

其中，表示第v+1步时三角形单元E_i的测度，表示第v步时三角形单元E_i的测度，即在新位置的单元测度等于旧单元测度与三条边扫过区域的测度之和。在单元新位置的守恒量平均值可更新为

$\left|E_i^{[v+1]}\right|U_i^{[v+1]}=\left|E_i^{\left[v\right]}\right|U_i^{\left[v\right]}+\underset{j=1}{\overset{3}{\mathrm{\&Sigma;}}}\hat{F}(U_{\mathrm{ij}}^L,U_{\mathrm{ij}}^R)---\left(22\right)$

这里通量采用如下的迎风格式

$\hat{F}(U_{\mathrm{ij}}^L,U_{\mathrm{ij}}^R)=\max \left\{\right|D_{\mathrm{ij}}|,0\}{U}_{\mathrm{ij}}^L+\min \left\{\right|D_{\mathrm{ij}}|,0\}{U}_{\mathrm{ij}}^R---\left(23\right)$

其中，和是由单元平均值利用SGM方法重构出来的边l_ij的左右极限值，SGM方法见式(9)。可见满足如下守恒性要求

$\hat{F}(U_{\mathrm{ij}}^L,U_{\mathrm{ij}}^R)=-\hat{F}(U_{\mathrm{ij}}^R,U_{\mathrm{ij}}^L)$

从而守恒性公式(19)成立，即在网格迭代第v步时的总质量与网格迭代第v+1步时的总质量相等。此外，格式能够精确捕捉静态。不妨令当网格迭代到第v步时，三角形单元Ei所在河床底部x、y方向的水流速度自由面高程const为一给定的常数。由SGM方法知初始数据重构出来的多项式(9)满足又由底部高程函数的线性离散知式(13)中的因此重构出来的单元边界处速度的左右极限都为0，代入式(23)知又由于式(21)可知，利用守恒型插值公式(22)更新之后的守恒量仍然满足 即自适应移动网格能保证静态流的捕捉。

为便于说明本发明效果起见，以下提供四个应用实例：

应用实例一：二维静态流问题

这个应用实例可以验证本发明所提供格式能够保证静态流。初始时刻水是静止的，即在x方向的速度u(x，y，0)＝0和在y方向的速度v(x，y，0)＝0。河床底部高程h(x，y)和初始时刻的水深B(x，y)用如下函数表示：

$\left\{\begin{array}{c}h(x,y)=1-B(x,y)\\ B(x,y)=0.8\exp \{-50[{(x-0.5)}^2+{(y-0.5)}^2\left]\right\},(x,y)\&Element;\left[\mathrm{0,1}\right]\&times;\left[\mathrm{0,1}\right]\end{array}\right.$

其中exp表示指数函数，(x，y)表示位置坐标。参数α＝50和β＝0.25。表1显示的是在t＝1.7分别用固定网格和和自适应移动网格下GRP格式的计算结果，结果表明，在非结构网格上GRP格式计算出来的自由面高程h+b，x、y方向的流速u、v与精确解的误差都处于机器误差范围内，所以该格式能够保证静态流。移动网格方法的计算结果同样也表明格式保证静态流。

表1  固定网格和和自适应移动网格下GRP格式的计算结果(t＝1.7)

应用实例二：平底二维溃坝问题

平底二维溃坝问题是非连续瞬间流的经典问题，可以用来检验本发明技术方案的激波扑捉能力。其中，河底底部高程h等于0，并假设在瞬间溃坝之前水是静止的，上游水库水位为10米，下游水位高度为5米。如图3a中箭头指出，区域中间坝的高度为75米，宽为10米。左右两边采用浸润边界条件，其他的边采用反射边界条件；计算区域为200×200米²。

在对区域进行初始网格生成两套完全独立的网格，程序在两套网格上分别运行，网格的单元数和节点数见表2，移动网格的是在粗网格的基础上进行移动。移动网格时的参数α＝20和β∈0.25。表2显示的是t＝7.2时粗网格固定网格1、密网格固定网格2和移动网格(AMM)的计算结果的网格和等值线，相应单元数，节点数，最长/最短边长、最大/最小单元面积和CPU计算时间如表2所示。

表2固定网格和自适应移动网格下平底二维溃坝数值计算的网格参数与计算时间

图3为本发明应用实例二中平地二维溃坝在t＝7.2时网格以及自由面高程h+b的等值线示意图，其中图3a为固定的粗网格，图3b为固定的密网格，图3c为移动的粗网格，图3d为在粗网格下的计算结果，图3e为在固定的密网格下的计算结果，图3f为在移动的粗网格下的计算结果。由图3可见，尽管粗固定网格(图3a)和自适应移动网格(图3c)的网格数一样，但两者在刻画高度H的梯度时却有很大差异。相比粗固定网格，自适应移动网格t＝7.2把有限的网格数集中在H梯度大的区域，包括两个非凸区域和右行激波，而在H梯度小的其他区域的网格相比固定网格却更为稀疏。正因为自适应移动网格的这个特性，在H梯度大的非凸区域和右行激波区域准确度和分辨率都相对更高(图3f)，在梯度小的其他区域却也能保证得到其准确度。

尽管细固定网格加密了4倍网格数14396，节点数7409，最小边长1.84m，最小面积1.71m²(表2)，但并没有因此而获得比自适应移动网格更高的准确度和分辨率图b，相反，细固定网格的CPU计算时间为67s，而自适应移动网格的CPU计算时间为40s，仅为前者的0.6。从其中一横断面y＝132米来看(图4)，更清晰的看到，在H变化大的区域自适应移动网格计算结果非常逼近细固定网格(图4)，相反，采用同等数量网格数的粗固定网格确在H梯度大的区域都与自适应移动网格和细固定网格有比较明显差异。图4是在直线y＝132上截取出来的速度图。其中实线表示的是在固定的密网格的计算结果，圆圈‘o’表示的是在移动的粗网格的计算结果，倒三角表示的是在固定的密网格上的计算结果。

此外，有意思的是，该算例存在2个非凸区域，但在整个时间尺度t＝0～7.2移动网格皆没有出现网格缠绕现象，而这对于实现移动网格在一般区域上求解浅水方程是非常重要的前提条件。

应用实例三：圆形溃坝问题

圆形溃坝问题可以验证本发明方法在激波扑捉方面的对称性。其中计算区域为50×50米²，河底底部高程h为0，圆形溃坝直径为20米。假设在瞬间溃坝的初始条件为坝内水位为2.5米坝外水位为0.5米。其边界条件全部采用浸润形式；其他参数见附图5。

在对区域进行初始网格生成两套完全独立的网格，程序在两套网格上分别运行，网格的单元数和节点数见表2，移动网格的是在粗网格的基础上进行移动。参数α＝40和β∈0.7。附图5显示的是t＝2时粗网格、密网格和移动网格的计算结果的网格和等值线，相应单元数，节点数，最长/最短边长、最大/最小单元面积和CPU计算时间如表3所示。

表3固定网格和自适应移动网格下圆形溃坝数值计算的网格参数与计算时间

与上一个溃坝例子一样，自适应移动网格求解圆形溃坝过程在梯度大的激波如r≈18米具有比粗固定网格更高的准确度和分辨率(图5)；相比细固定网格网格数23198，节点数11800，最小边长0.32米，最小面积0.05米²(表3)，自适应移动网格计算的H具有相近的准确度和分辨率；但从计算成本来看，密固定网格的CPU计算时间为61s，而自适应移动网格的计算时间为38s，也仅为前者的0.6。从其中一横断面y＝25米来看，更清晰的看到，自适应移动网格计算结果在H梯度大的区域非常逼近细固定网格(图6)；然而，采用粗固定网格确很难扑捉激波区域。图6是在直线y＝25上截取出来的自由面高程图。其中实线表示的是在固定的密网格的计算结20111107果，圆圈‘o’表示的是在移动的粗网格的计算结果，倒三角表示的是在固定的密网格上的计算结果。

更为重要的是，不管粗固定网格还是加密4倍的密固定网格，自适应移动网格都相对更好的表征圆形溃坝过程的对称性现象。甚至还尝试了加密固定网格16倍(表3)，也并没有因为网格数的增加而有效改善对称性激波扑捉能力。为了更为清晰说明移动网格描述上述对称性现象的能力，本发明列举出来t＝0，0.4，0.8，1.2，1.6(t＝2.0见图5)的溃坝过程及相应的移动网格。图7为本发明应用实例三中圆形维溃坝在t＝0，0.4，0.8，1.2，1.6时的网格变化示意图，其中图7a为t＝0时的网格图，图7b是t＝0时相应的自由面高程h+b的等值线示意图，图7c为t＝0.4时的网格图，图7d是t＝0.4时相应的自由面高程h+b的等值线示意图，图7e为t＝0.8时的网格图，图7f是t＝0.8时相应的自由面高程h+b的等值线示意图，图7g为t＝1.2时的网格图，图7h是t＝1.2时相应的自由面高程h+b的等值线示意图，图7i为t＝1.6时的网格图，图7j是t＝1.6时相应的自由面高程h+b的等值线示意图。由图7可知，自适应移动网格是动态地刻画或逼近物理量H的梯度分布和激波产生、移动过程。从计算原理来看，主要是因为自适应移动网格不仅依据变分原理在激波区域加密了网格，同时能使得网格的空间分布与H值的空间分布具有平行特征；不仅如此，在网格的移动过程中的守恒量在新网格上的插值采用的是基于几何的方法，近似于一种积分平均，此守恒性插值对解具有磨光效应。

Claims (5)

1.一种复杂边界及实际地形上溃坝洪水的自适应模拟方法，其特征在于：将待模拟区域按网格进行处理，每个网格是一个三角形单元，通过守恒型的物理量U反映每个网格上水量信息，物理量U的向量形式如下，

$U=\left[\begin{array}{c}h+b\\ \mathrm{hu}\\ \mathrm{hv}\end{array}\right]$

其中，h表示水深，b表示河床底部高程，u表示坐标轴x方向的流速，v表示坐标轴y方向的流速；

按网格处理的具体实现包括以下步骤，

步骤1，根据待模拟区域，输入物理区域的初始网格和逻辑区域的初始网格；在时刻t_n，物理区域中第i个三角形单元E_i所在网格记为物理区域中第i个三角形单元E_i上物理量U的积分平均值记为令n=0，计算在初始时刻t_n＝t₀的积分平均值

步骤2，进行移动网格，调整物理区域上的网格位置，移动任一网格的具体实现包括以下步骤，

步骤2.1，令迭代步数ν＝0，ν＝0时网格等于当前时刻t_n的网格ν＝0时积分平均值等于当前时刻t_n的积分平均值

步骤2.2，计算表征物理量变化梯度大小的控制函数ω_i，计算公式如下，

$\left\{\begin{array}{c}{\mathrm{\&omega;}}_i=\sqrt{1+\mathrm{\&alpha;}{\widetilde{\mathrm{\&omega;}}}_i^2(\mathrm{\&beta;},h+b)}\\ {\widetilde{\mathrm{\&omega;}}}_i(\mathrm{\&beta;},h+b)=\min \{1,|{\&dtri;}_{\mathrm{\&xi;}}(h+b){|}_i/\mathrm{\&beta;}\underset{i}{\max }\left|{\&dtri;}_{\mathrm{\&xi;}}(h+b){|}_i\right\}\end{array}\right.$

其中算子▽_ξ表示在参考坐标下的梯度，参考坐标为逻辑区域的坐标；h表示水深，b表示河床底部高程，α和β为正参数；

步骤2.3，根据步骤2.2所得控制函数ω_i，通过以逻辑区域作为中间变量，移动网格至新的位置，新的位置记为网格

步骤2.4，在网格上更新迭代到第v+1步时的积分平均值包括通过守恒型插值公式得到积分平均值守恒型插值公式如下，

$\left|E_i^{[v+1]}\right|U_i^{[v+1]}=\left|E_i^{\left[v\right]}\right|U_i^{\left[v\right]}+\underset{j=1}{\overset{3}{\mathrm{\&Sigma;}}}\hat{F}(U_{\mathrm{ij}}^L,U_{\mathrm{ij}}^R)$

其中，表示第v+1步时三角形单元E_i的测度，表示第v步时三角形单元E_i的测度，通量采用如下的迎风格式

$\hat{F}(U_{\mathrm{ij}}^L,U_{\mathrm{ij}}^R)=\max \left\{\right|D_{\mathrm{ij}}|,0\}{U}_{\mathrm{ij}}^L+\min \left\{\right|D_{\mathrm{ij}}|,0\}{U}_{\mathrm{ij}}^R$

其中，和是由单元平均值利用SGM方法重构出来的边l_ij的左右极限值，|D_ij|是网格移动到网格时边l_ij所扫过区域D_ij的测度，边l_ij表示三角形单元E_i的边；

步骤2.5，判断是否满足预设的迭代结束条件，如果是则迭代结束，更新在下一时刻t_n+1的网格位置到网格更新当前时刻t_n的积分平均值转入步骤3；否则令迭代步数ν＝ν+1，转步骤2.3；

步骤3，根据步骤2的移动网格结果，求取下一时刻t_n+1的网格和积分平均值包括以下步骤，

步骤3.1，在网格上，以步骤2更新的积分平均值作为初始值，求取下一时刻t_n+1的积分平均值的近似解；

步骤3.2，如果下一时刻t_n+1＜TSTOP，TSTOP为预设的停止时刻，则令n＝n+1，以下一时刻为当前时刻返回执行步骤2，否则输出下一时刻t_n+1的网格所在位置和积分平均值作为模拟结果。

2.根据权利要求1所述复杂边界及实际地形上溃坝洪水的自适应模拟方法，其特征在于：步骤2.5中，采用或迭代步数ν＝5作为预设的迭代结束条件。

3.根据权利要求1或2所述复杂边界及实际地形上溃坝洪水的自适应模拟方法，其特征在于：步骤2.3中，通过解离散的网格方程移动网格至新的位置，离散的网格方程如下，

其中，表示迭代的网格坐标位置的中间变量，表示迭代时的权重，μ_i为松弛因子，表示当网格x_i迭代到第ν步时的邻居节点，N(i)表示网格x_i周围邻居节点的个数，且

其中，|ξ_ijξ_i|,|l_ij|分别表示逻辑区域上的网格边长及对偶网格边长，Δτ和σ是控制网格质量的变量，ω_ij是边l_ij的两个相邻单元上的控制函数w_i与w_j的平均值；|V_i|表示逻辑区域上以节点ξ_i为中心的对偶单元V_i的面积。

4.根据权利要求1或2所述复杂边界及实际地形上溃坝洪水的自适应模拟方法，其特征在于：

步骤3.1中，求取积分平均值的实现方式如下，

首先，以步骤2所得积分平均值作为初始值，利用自由面梯度法SGM重构分片线性的初始函数，初始函数见下式

其中，W(x,y)表示向量W的线性逼近函数，是向量W于当前时刻t_n在三角形单元E_i上的近似值，向量W＝(h+b,u,v)^T；r=(x-x₀,y-y₀)表示位置坐标(x,y)相对三角形单元E_i中心(x₀,y₀)的相对坐标，梯度是由三角形单元E_i上的三个相邻单元E_i1、E_i2、E_i3上向量W的值，是一个选定的限制器；

然后，将在中点时刻边l_ij的中点处物理量代入下式，

$S{\left(U\right)}_i^{n+\frac{1}{2}}={[0,{\mathrm{gb}}_i{(h_x+b_x)}_i^{n+\frac{1}{2}},{\mathrm{gb}}_i{(h_y+b_y)}_i^{n+\frac{1}{2}}]}^T$

得到物理量U的源项S(U)的离散， $S\left(U\right)=\begin{array}{}\end{array}\left[\begin{array}{c}0\\ \mathrm{gb}{(h+b)}_x\\ \mathrm{gb}{(h+b)}_y\end{array}\right],$ ((h+b)_x,(h+b)_y)表示自由面高程的梯度，

其中，表示源项S(U)的近似值，g表示重力加速度，b_i表示在三角形单元E_i中心处河床底部高程；上标表示中点时刻Δt＝t_n+1-t_n表示时间步长；(b_x,b_y)表示河床底部高程b的梯度；h_x表示水深h的x方向变化率，h_y表示水深h在y方向的变化率；

最后，利用下式得到下一时刻t_n+1的积分平均值的近似解，

$U_i^{n+1}=U_i^n-\frac{\mathrm{\&Delta;t}}{\left|E_i\right|}\underset{j=1}{\overset{3}{\mathrm{\&Sigma;}}}{F}_{n_{\mathrm{ij}}}{\left(U\right)}_{\mathrm{ij}}^{n+\frac{1}{2}}{\mathrm{\&Delta;l}}_{\mathrm{ij}}+\mathrm{\&Delta;tS}{\left(U\right)}_i^{n+\frac{1}{2}}$

其中，Δt＝t_n+1-t_n表示时间步长，|E_i|表示三角形单元E_i的面积，Δl_ij为边l_ij的长度；和分别表示通量F(U)和与源项S(U)的近似值，通量 $F\left(U\right)=\left[\begin{array}{c}\mathrm{hu}\\ {\mathrm{hu}}^2+\frac{1}{2}{g(h+b)}^2\\ \mathrm{huv}\end{array}\right].$

5.根据权利要求4所述复杂边界及实际地形上溃坝洪水的自适应模拟方法，其特征在于：在中点时刻边l_ij的中点处的物理量通过解广义黎曼问题得到近似值，具体方式如下，

$U_{\mathrm{ij}}^{n+\frac{1}{2}}=U_{\mathrm{ij}}^n+\frac{\mathrm{\&Delta;t}}{2}{\left(\frac{\&PartialD;U}{\&PartialD;t}\right)}_{\mathrm{ij}}^n$

在当前时刻t_n边l_ij的中点处的物理量通过下式求出近似值

$U_{\mathrm{ij}}^n=\underset{t\&RightArrow;0^{+}}{\lim }U(0,t)=\mathrm{RM}\_\mathrm{Solver}(0;U_L,U_R)$

其中，RM_Solver(0;U_L,U_R)表示以U_L和U_R作为左右状态的黎曼问题解在时间轴上的值，U_L和U_R表示由SGM方法重构出来在单元边界的左右极限值；

时间方向导数通过以下关于水深时间变化率和水速时间变化率的二元一次方程组求出，

$\left\{\begin{array}{c}{a}_{-}^1{\left(\frac{\&PartialD;h}{\&PartialD;t}\right)}_{\mathrm{ij}}^n+a_{-}^2{\left(\frac{\&PartialD;u}{\&PartialD;t}\right)}_{\mathrm{ij}}^n=d_{-},\\ a_{+}^1{\left(\frac{\&PartialD;h}{\&PartialD;t}\right)}_{\mathrm{ij}}^n+a_{+}^2{\left(\frac{\&PartialD;u}{\&PartialD;t}\right)}_{\mathrm{ij}}^n=d_{+},\end{array}\right.$

其中，和是常数。