2.2　数值计算方法

采用边界拟合能力强和易于局部网格加密的三角形和四边形网格剖分计算域，利用基于水位-体积关系的斜底单元模型，有效解决了小尺度线状地形模拟难题；以能够有效捕获激波的Godunov型有限体积法为框架，运用在时间上和空间上均具有二阶精度的MUSCL-Hancock预测-校正格式离散洪水控制方程，采用HLLC近似Riemann算子计算对流数值通量，采用直接近似方法计算扩散数值通量，并结合斜率限制器以保证模型的高分辨率特性，避免在间断或大梯度解附近产生非物理虚假振荡；基于单元中心型底坡项近似方法，在不使用任何额外动量通量校正项的前提下模型能保持通量梯度与底坡项之间的平衡，即模型具有和谐性质；采用半隐式格式处理摩阻项，该半隐式格式既能保证不改变流速分量的方向，也能避免小水深引起的非物理大流速问题，有利于计算稳定性；实现了固壁、水位、流量、自由出流等边界条件；基于克朗稳定条件实现了数值模型的自适应时间步长技术。

2.2.1　计算网格

鉴于非结构三角形、四边形混合单元具有复杂边界拟合能力强、便于网格生成和局部加密等特点，Hydro MPM2D_FLOW采用非结构三角形、四边形单元作为计算网格。以三角形单元为例（四边形网格类似），网格拓扑结构如图2.3所示。其中Ci为待计算单元，其顶点1—顶点2—顶点3排序服从逆时针方向；与顶点k相对的边为Γi，k，其外法向单位向量为ni，k；单元Ci的邻接单元中，与顶点k相对的单元为Ci，k。

对于每一个节点而言，其网格拓扑信息包括节点序号，节点坐标，节点周围的单元；对于每一条边而言，其网格拓扑信息包括边序号，边的左、右单元，边的始、末节点；对于每一个单元而言，其网格拓扑信息包括单元序号，单元的三个顶点，单元的三条边，单元的三个邻接单元。在模型开始计算之前，需要根据网格文件构造包含上述拓扑信息的计算网格系统。

2.2.2　斜底单元模型

在计算网格结构中，有两种底高程定义方式：①将底高程定义于单元形心处，单元内的底高程为均一值；②将底高程定义于单元顶点处，单元内的底高程服从线性分布。在地形表达精度方面，第一种方式仅为一阶精度，而第二种方式具有二阶精度。此外，在水动力数值模拟的实际工程中，往往包含由生产堤和公路等线状建筑物组成的奇异地形。该类奇异地形具有低水位干出、高水位淹没的性质，必须在模型中予以准确表达。然而，如果采取第一种底高程定义方式，由于此类奇异地形的空间尺度要远小于满足计算效率要求的网格尺度，因此需要采用局部网格加密方法以表达此类奇异地形，此时不仅网格数量剧增，而且由于稳定条件的限制，小尺度网格将导致模型的计算时间步长大幅度减小，严重影响模拟效率。另一方面，如果采取第二种底高程定义方式，在网格划分之前，将一系列节点预先布置在生产堤和公路等线状建筑物上，进而使该类奇异地形在网格系统中以“边”的形式得以表达。高水位时，奇异地形所在网格边被淹没过水；低水位时，奇异地形所在网格边的物质通量为零，起到阻水作用，故第二种底高程定义方式可以实现该类奇异地形的准确模拟。综合上述两方面的原因，Hydro MPM2D_FLOW采用第二种底高程定义方式，即斜底单元模型。

以三角单元为例（四边形单元类似），在斜底单元模型中，守恒变量h代表单元平均水深，单元的水量为hΩ，其中Ω为单元面积；水位η代表单元内含水部分的水面高程，且假设单元内含水部分的水面为一平面（图2.4）。其中，斜线阴影面为水面，三角形123为单元底面，图2.4（a）中三个顶点的水深均大于零，图2.4（b）和图2.4（c）存在水深为零的顶点。

不失一般性，假设单元Ci三个顶点的底高程bi，1、bi，2和bi，3满足关系：bi，1≤bi，2≤bi，3。由水量守恒原理，可得单元水深与单元水位之间的转换计算式。

1.已知水深，计算水位

其中　　γ1=bi，3-3bi，1

2.已知水位，计算水深

若将底高程定义于单元形心处，且假设单元内的底高程为均一值，则单元水位的最小值为单元形心处的底高程值，在干湿界面计算时需要重构干单元底高程。然而，由式（2.17）可知，斜底单元模型的引入，使单元水位的最小值降低为单元三个顶点底高程的最小值，避免了干湿界面计算时重构干单元底高程，提高了模型干湿界面处理能力，且有利于设计具有和谐性的计算格式。

在斜底单元模型中，单元具有全淹没、局部淹没和全干三种状态。以三角形单元为例（四边形单元类似），Hydro MPM2D_FLOW采用如下状态判别方法：

对于三角形单元，由于底高程定义于单元顶点处，根据不共线三点决定一个平面的原理可计算出单元底高程的平面方程，继而可得到该平面在x和y方向的斜率，即单元的底坡斜率：

式中：xi，k、yi，k、bi，k分别为单元Ci第k个顶点的x方向坐标、y方向坐标和底高程（k=1，2，3）；Ωi为单元Ci的面积。

对于四边形单元，单元的底坡斜率可根据单元内2个三角形底坡斜率的面积加权平均得到。

2.2.3　有限体积离散

在任意控制体Ω上对式（2.1）所示的控制方程进行积分得

对式（2.22）运用Green公式将体积分转化为沿其边界的线积分，可得

式中：∂Ω为控制体Ω的边界；n为边界∂Ω的外法向单位向量；dΩ、dl分别为面积微元和弧微元；Fadv=［Eadv，Gadv］T；Fdiff=［Ediff，Gdiff］T；Fdis=［Edis，Gdis］T。

在平面二维网格中，线积分可分别由下式计算：

式中：N为网格的边数（三角形单元N=3，四边形单元N=4）。假设Ui为U在单元Ci内的平均值，即

则有

式中：Ωi为单元Ci的面积；、、、ni，k、Li，k分别代表单元Ci第k条边的对流数值通量、雷诺应力引起的扩散数值通量、二次流引起的扩散数值通量、外法向单位向量和长度；Si为源项近似。

2.2.4　数值通量计算

由于积分平均，物理变量在每个单元内部为常数，在整个计算域内形成阶梯状分布，因此在单元界面处物理量存在间断，即界面左、右两侧的物理量不相等，故而在界面处构成了一个局部Riemann问题。通过Riemann问题的求解可得到界面处的对流数值通量，即

式中：UL、UR分别为界面左侧和右侧的守恒向量。

由浅水方程的旋转不变性可得：

式中：T和T-1分别为旋转矩阵及其逆矩阵。

定义：

式中：u⊥、u∥分别为与界面垂直和平行的流速分量。

则有

由式（2.34）可知，二维浅水方程的对流数值通量计算可转化为界面处一维Riemann问题求解。

如图2.5所示，在ξ-t平面上，二维浅水方程Riemann解中存在三个波，分为三类：激波（shock wave）、稀疏波（rarefaction wave）和接触波（contact wave）。三个波将初始的两个常状态区分割为四个常状态区，其中左、右两个区域的状态分别为初始的左、右状态值，中间两个区域的状态则由Riemann解的结构确定。

目前较常用的近似Riemann求解器主要有FVS格式、FDS格式、Osher格式、Roe格式、HLL格式、HLLC格式等。由于HLLC格式满足熵条件，且在合理计算波速的情况下适应干湿界面计算，因此Hydro MPM2D_FLOW采用该格式计算二维浅水方程的对流数值通量。

式中：，均由式（2.34）计算；UL、UR分别为局部Riemann问题所在界面左侧和右侧的守恒向量；分别为Riemann解中间区域接触波左、右侧的数值通量；s1、s2、s3分别为左波、接触波和右波的波速（图2.5）。

接触波左、右侧的数值通量分别由式（2.36）和式（2.37）计算：

式中：分别为运用HLL格式计算得到的法向数值通量的第一、第二个分量。

运用HLLC求解器计算数值通量的关键在于波速近似。Hydro MPM2D_FLOW采用Einfeldt波速计算式：

式中：h*和u⊥，*为Roe平均。

由于激波的波速小于激波后面区域的特征速度，此时Einfeldt波速即为激波波速的Roe近似，因此使用Einfeldt波速可获得在激波附近更加准确的数值解。

由式（2.43）计算接触波的波速：

式（2.35）～式（2.43）即为对流数值通量的计算方法。

雷诺应力引起的扩散数值通量计算过程较为简单，即

式中：hi，k为单元Ci第k条边所在界面处的水深；为流速分量在单元Ci第k条边所在界面处的斜率。

在计算扩散数值通量时，界面处的水深值取界面两侧单元水深重构值的平均：

式中：、分别为单元Ci第k条边所在界面左、右侧的水深重构值。

界面处的流速梯度取界面两侧单元流速梯度的面积加权平均：

式中：、分别为单元Ci第k条边所在界面左、右侧单元的面积；、为界面左侧单元的流速分量斜率；分别为界面右侧单元的流速分量斜率。

二次流引起的扩散通量按式（2.48）计算：

式中：Di，k为界面处的二次流扩散应力项，取界面两侧单元的面积加权平均。

2.2.5　变量重构和斜率限制函数

如果物理量在单元内近似为常数分布，并且以各单元形心处的值作为界面处局部Riemann问题的初始间断条件，则相应的计算格式在空间上仅具有一阶精度，存在较大的数值耗散。对水流模拟而言，一阶精度的计算格式能基本满足工程应用的要求。然而，对泥沙运动或者污染物扩散模拟而言，一阶格式的过度数值耗散往往导致计算结果失真。为了提高格式的空间精度，在构造界面处局部Riemann问题时，需要采用比分片常数近似函数更高精度的重构方法对界面左、右两侧的变量进行重构，并基于重构变量求解界面处的数值通量。Hydro MPM2D_FLOW针对不同的单元干湿状态，采用不同的变量重构方法。

若界面某侧的单元为局部淹没状态，或者该侧单元存在处于局部淹没状态的邻接单元，则界面该侧的流速重构值为该侧单元形心处的值：

式中：为界面重构值；ui、vi为单元形心处的值。

同时，界面该侧的水深重构值按式（2.51）计算：

式中：为界面重构值；ηi为单元的水位值；分别为单元Ci第k条边首尾两个节点底高程的最小值和最大值。

若界面某侧的单元为全淹没状态，且该侧单元的所有邻接单元均为全淹没状态，则通过线性重构并结合限制函数计算界面该侧的水位、流速重构值：

式中：p为重构变量，如η、u或v；（x，y）为界面中点处的坐标；r为界面中点相对于单元形心的位置矢量；为限制梯度。

式中：为原始梯度；φi为限制函数。

式中：为单元顶点处未受限制的重构值，ri，k为单元顶点相对于单元形心的位置矢量；和分别为当前单元及其邻接单元形心处变量的最大值和最小值，即；pi为单元Ci形心处的变量值。

对于三角形单元，假设pi，1、pi，2、pi，3分别为变量p在单元Ci三个邻接单元形心（xi，1，yi，1）、（xi，2，yi，2）和（xi，3，yi，3）处的值，则变量p在单元Ci的梯度计算式为

对于四边形单元，变量p在单元Ci的梯度由四边形内2个三角形单元梯度的面积加权得到。

水深重构值等于水位重构值减去底高程：

2.2.6　底坡项近似

在每一个计算步，Hydro MPM2D_FLOW均对所有单元进行连续方程和动量方程的更新，即同时更新单元的水深和流速。因此，需要对所有单元进行底坡项处理。

若单元为全淹没状态，Hydro MPM2D_FLOW采用单元中心型近似方法处理底坡项：

否则，Hydro MPM2D_FLOW采用DFB（Divergence Form for Bed slope source term）技术处理底坡项：

式中：N为单元边数；bi，k为单元Ci第k条边的中点位置底高程。

采用上述底坡项近似方法，保证了模型的和谐性，证明过程如下。

考虑静水条件，即流速为零、水位为常数。由于界面两侧的水位相等、流速为零，所以所有界面的物质通量为零。单元Ci所有边的外法向对流数值通量为

式中：i、k分别为单元序号和边序号；；h、z分别为边中点处的水深重构值和底高程。

若单元Ci为全淹没状态，则

式中：η为静水条件下的水位统一值。

将式（2.62）代入式（2.61），可得

以三角形单元为例（四边形单元类似），假设单元Ci第k个顶点的坐标及底高程分别为xi，j、yi，j和bi，j，顶点1—顶点2-顶点3服从逆时针排序，与顶点k相对的边为单元Ci第k条边，则有如下关系成立：

故有

因此，单元Ci三条边的外法向对流数值通量可化简为

同时，由于静水条件下流速为零，故扩散数值通量和摩阻项为零。因此有

若单元为半淹没单元，则由式（2.59）～式（2.61）可知，单元Ci三条边的外法向对流数值通量之和与底坡项近似值相等，因此，式（2.69）仍然成立。

综上所述，模型能完全保持水流的静止状态，即在不使用任何额外校正项的前提下，Hydro MPM2D_FLOW的二维水流模型具有和谐性。

2.2.7　摩阻项处理

一般地，摩阻项通过算子分裂法进行处理：

由于在摩阻项处理过程中，单元水深保持不变，即，因此，式（2.70）可简化为

复杂地形的陡峭坡面，使局部区域的水深较小、流速较大，使得式（2.71）对应的常微分方程系统的Lipschitz常数很大，因此，摩阻项可能引起刚性问题。此时，若采用一般的显式数值方法处理摩阻项，将显著影响数值计算的稳定性，或将极大减小全局时间步长，从而严重降低计算效率。为解决上述问题，需要采用隐式或半隐式格式处理摩阻项。然而，由于水深变量位于摩阻项的分母，一般的隐式或半隐式计算格式仍面临一些问题，如产生错误的大流速、改变流速分量的方向等。

综合考虑摩阻项处理的稳定性和计算效率，Hydro MPM2D_FLOW采用如下半隐式格式处理摩阻项。

定义：

则由式（2.71）可得

利用半隐式格式求解式（2.73），可得

即

同理可得

式中为利用数值通量对n时刻已知量进行更新得到的值。由式（2.75）和式（2.76）可知，采用的半隐式格式能保证不改变流速分量的方向，有利于计算稳定性。

2.2.8　柯氏力项处理

柯氏力项通过算子分裂法进行处理：

由于在柯氏力项处理过程中，单元水深保持不变，即，因此式（2.77）可简化为

取W=u+vi，将式（2.78）乘以虚数单位i，可得

对于欧拉前差，则有

误差方程为

对复数取绝对值，则有

因此，对于柯氏力项，欧拉前差是无条件不稳定的。同理可知，欧拉后差是无条件稳定的。因此，Hydro MPM2D_FLOW采用欧拉后差格式处理柯氏力项：

由式（2.82）可得

式中：和为考虑柯氏力项之前的流速值。

2.2.9　风应力和波浪辐射应力处理

采用直接近似法和显式格式处理风应力和波浪辐射应力：

式中：风应力和波浪辐射应力均为n时刻的单元中心近似值。

2.2.10　时间二阶积分

在二维浅水方程数值求解过程中，时间上离散可采用Runge-Kutta、Hancock预测-校正等多步格式，以提高模型的时间精度。从计算效率看，在一个时间步长内，两步Runge-Kutta法分别对所有单元界面计算两次黎曼问题，而Hancock预测-校正法只需要计算一次黎曼问题，因此，Hancock预测-校正法的计算效率相对较高；从格式稳定性看，Runge-Kutta法可满足TVD性质，而MUSCL-Hancock法在合理选择时间步长的前提下满足L1稳定。考虑模型的计算效率，Hydro MPM2D_FLOW采用Hancock预测-校正格式实现二维浅水方程数值求解的时间二阶积分。

2.2.10.1　预测步

在预测步，计算域的水流状态由t时刻更新至t+Δt/2时刻，其中Δt为计算时间步长。Δt可以设定为常数，也可以根据稳定条件进行自适应调整。预测步不考虑紊动扩散项。

若单元为全干或局部淹没状态，则以当前时刻的水位、流速作为预测步的结果：

式中：分别为t时刻的水位、流速；为预测步结果。

若单元为全淹没状态，则以水位、流速作为预测变量：

由于静水条件下流速为零、水位为常数，因此，由式（2.86）～式（2.89）可知，预测步的水流静止状态得以维持。

2.2.10.2　校正步

在校正步，计算域的水流状态由t时刻更新至t+Δt时刻。基于守恒形式浅水方程，由式（2.90）对单元值进行更新：

式中分别为单元Ci在t和t+Δt时刻的单元平均守恒向量；、分别为基于预测步结果的对流数值通量和扩散数值通量；Si为源项近似。

2.2.11　边界条件

一般情况下，模型的边界条件实现方式有两种：镜像单元法和直接计算数值通量法。其中，前者在基于结构网格的数学模型中应用较广，而后者被广泛运用于基于非结构网格的数值模型。Hydro MPM2D_FLOW采用直接计算数值通量的方式实现边界条件。

边界处的对流数值通量为

式中：h*、u*、u⊥，*分别代表边界边中点处的水深、x方向流速分量、界面外法向方向流速分量；b为边界边中点处的底高程值；nx、ny分别为边界边的外法向单位向量在x和y方向的分量。

由式（2.91）可知，计算边界边对流数值通量的关键在于估算边界边中点处的水深、流速等物理量。

2.2.11.1　急流开边界条件

由浅水方程特征理论可知，水流为急流状态时信息沿下游传播，因此急流出口边界的扰动不会对计算域内的水流状态产生影响，故急流边界处的水深、流速等水力要素值采用边界边所在内单元的水力要素值：

在数值模拟中，常用自由出流边界条件，其实现方式与急流开边界条件相同。

2.2.11.2　缓流开边界条件

假设边界边的左侧位于计算域之内，右侧位于计算域之外，基于浅水方程特征理论，由一维浅水方程的Riemann不变量可得：

式中：hL、u⊥，L分别为边界边所在内单元形心处的水深和界面外法向方向流速分量。

式（2.94）中存在两个未知变量，因此需要结合边界条件建立关于h*和u⊥，*的等式，进而通过求解方程组计算h*和u⊥，*。

（1）水位边界条件。在边界上给定水位η=η（t），则边界上水深为h*=η（t）-b。由式（2.94）可得

（2）流速边界条件。在边界上给定边界边外法向方向的流速分量u⊥，*，由式（2.95）可得

（3）单宽流量边界条件。在边界上给定边界边的外法向方向单宽流量q⊥，*=q（t）：

由式（2.97）有

将式（2.98）代入式（2.94）可得

式中：

当q⊥，*＜0时，为入流边界条件。通过对函数，*进行分析可知，如图2.6所示，；当aL＜0时，f（c*）为单调递增函数；当aL＞0时，在区间（0，aL/3］内为单调递减函数，在区间［aL/3，+∞］内为单调递增函数。因此，入流边界条件时，方程f（c*）=0存在唯一的实数根。

当q⊥，*＞0时，为出流边界条件。通过对函数，*进行分析可知，如图2.7所示，；当aL＜0时，f（c*）为单调递增函数，f（c*）=0无实数根；当0＜aL＜3（gq⊥，*）1/3时，f（c*）的最小值大于零，f（c*）=0无实数根；当aL=3（gq⊥，*）1/3时，f（c*）的最小值等于零，f（c*）=0有且仅有一个实数根；当aL＞3（gq⊥，*）1/3时，f（c*）=0有两个互不相等的实数根。

值得注意的是，出流边界条件时，图2.7所示四种情况并非全部成立。由于出流边界处于缓流状态，即

定义：

则φ（h*）的导数为

则有

故可得h*的取值范围为

将代入式（2.101）可得

因此

即

故有

因此，出流边界条件时，图2.7所示四种情况中，仅有aL＞3（gq⊥，*）1/3情况成立。虽然此时f（c*）=0有两个互不相等的实数根，但仅有一个实数根满足式（2.100）所示缓流条件。

采用直接求根的方法求解式（2.99）所示的一元三次方程，并结合上述可行解范围分析，对求取的方程根进行取舍。

（4）断面流量边界条件。当模拟实际河道上的洪水演进时，往往给定上游断面的流量过程。由于上游断面被淹部分较宽，水流存在一定的横比降，因此不能简单地通过将上游断面的流量除以断面被淹部分宽度实现单宽流量计算。

Hydro MPM2D_FLOW根据边界单元的过水断面面积正比例分配断面总流量：

式中：l为边界单元的边长度；Qi为边界断面的第i条边的流量。

第i条边的单宽流量为

通过式（2.110）计算得到各条边的单宽流量，进而结合单宽流量边界条件实现方法，进行相应的边界条件计算。

为获取边界处的x和y方向流速分量，还需要计算边界处的切向流速。假设边界处的切向流速等于边界边所在内单元形心处的切向流速，即

通过式（2.112）将局部坐标系下的流速值转换为笛卡尔坐标系下的流速值：

基于h*、u⊥，*、u*和v*计算结果，通过式（2.91）计算边界边的对流数值通量。

2.2.11.3　固壁边界条件

在固壁边界处，边界上的法向流速为零，且假设边界上的水深值等于边界边所在内单元形心处的水深值，即

结合式（2.91），边界边的对流数值通量为

2.2.12　稳定条件

由于采用显式格式求解浅水方程，为保持格式的稳定，时间步长Δt受克朗稳定条件的限制：

式中：Δt为时间步长；Cr为克朗（Courant）数，0＜Cr≤1，一般情况下取Cr=0.8；u⊥和h为界面的Roe平均；N为计算网格的单元总数。

2.2.13　数值计算流程

二维浅水方程数值计算流程如图2.8所示。