三角形网格下二维浅水方程的高分辨率格式

A Non-Oscillatory Scheme for Shallow Water Equations on Triangular Meshes

摘要: 对于二维浅水方程问题，本文在三角网格下，基于CBC (Convection Boundedness Criterion)准则采用有限体积法离散建立了新的高分辨率NVSF (Normalized Variable and Space Formulation)格式，NVSF格式是为解决NVF高分辨率格式在不规则区域的应用所提出的，文章通过一些典型算例的精确解与数值解的比较，表明新构造的数值格式具有二阶精度，并且与准确解有很好的逼近效果，能很好地抑制在间断解处的非物理震荡。

Abstract: For the shallow water equations, a new finite volume NVSF (Normalized Variable and Space For-mulation) scheme is constructed under the triangular mesh. This scheme is based on the CBC (Convection Boundedness Criterion) criterion. NVSF scheme is adopted to solve the application of NVF high resolution scheme on triangular meshes, by comparing the exact solution of the typical example with the numerical solution. It is shown that the new numerical scheme has second order accuracy and has a good approximation near the exact solution. It can suppress the unphysical os-cillation of the discontinuity.

1. 引言

浅水方程数值方法的研究一直是计算流体力学研究的主要方向。对于间断问题的求解是一个比较难解决的问题，当前求解二维浅水方程的数值计算方法有很多，如有限差分法、有限元法、特征法、有限体积法等。有限差分法在处理问题时相对比较简单效率较高，但是在处理复杂区域问题时较为困难；而有限元在处理复杂边界问题时具有很大的优势，但是在处理大梯度问题时是比较困难的，有限体积法是基于无结构网格下，处理强间断问题或大梯度问题的一种高效的方法。因此，目前求解二维浅水方程的最广泛的方法是有限体积法。而且在模拟浅水方程问题上已经取得了一些成果。如前面已经有学者基于特征思想的高阶思想来求解浅水方程 [1] ；王立辉等提出的Lax-Wendroff格式和Lax-Friedrichs格式交替使用形成的复合型有限体积格式 [2] ；窦红、汪继文提出的将一种Van Albada型可微的限制器函数引入二维浅水方程 [3] ；王昆等提出了一种基于HLL格式的近似Riemann解为基础的重构的Godunov格式 [4] ；朱华君、宋松和提出的将非结构ENO型有限体积法运用到二维浅水方程的求解中 [5] 。这些格式都是基于TVD方法来求解二维浅水方程的，而且对于浅水方程的求解大部分都是基于TVD格式的。本文中我们基于CBC准则来构造一种新的高分辨率格式来求解二维浅水方程。

对于间断问题的求解是一个比较难解决的问题，在一致网格下，人们已经构造了许多经典的数值格式，如一阶迎风格式FOU (First-Order Upwind)、Lax-Friedrichs [6] 格式等低阶格式，这种低阶格式虽然能使计算稳定且不产生震荡，但是数值解可能会在间断处存在数值耗散，因此又构造了一些其他的重要格式，如CD (Central Difference)、SOU (Second-order Upwind)等高阶数值格式，虽然这些新构造的数值格式具有高阶精度，但是这些数值格式得到的数值解可能会产生虚假震荡导致图像失真。将这些数值格式和对流有界性准则相结合构造出NVF (Normalized Variable Formulation) [7] 高分辨率格式，如HOAB [8] ，SMART [9] ，OSHER [10] 等高阶格式。但是，许多二维问题具有不规则边界，NVF高分辨率格式无法直接使用，因此将NVF高分辨率格式推广到非结构网格下对于数值计算十分重要。M. S. Darwish [11] 提出的非结构网格下NVSF (Normalized Variable and Space Formulation)格式是这方面的基础研究工作之一。由于CBC格式能保证数值格式的有界性，本文我们提出一种新的高分辨率格式来求解二维浅水方程问题的对流项。基于三角形网格下，在M. S. Darwish [11] 的研究基础上，结合CBC (Convection Boundedness Criterion) [12] 对流有界性准则构造一个新的NVSF (Normalized Variable and Space Formulation)格式，通过数值算例，给出数值解与准确解的比较，可以得到该格式对于准确解具有良好的逼近效果。本文共有六节组成，具体安排如下：第一节为引言部分，第二节为非结构网格下新格式的构造，第三节为二维浅水方程，第四节为时间离散方法，第五节为给定具体的数值算例进行分析，第六节给出结论。

2. 非结构网格下新格式的构造

2.1. 非一致网格下NVSF的正则化

在非一致网格下，界面f的值不仅与相邻的节点值 ${\varphi }_U,{\varphi }_C,{\varphi }_D$ 有关，而且也与这些节点所在的位置有关，如图1所示，从而可以得到 ${\varphi }_f=f\left({\varphi }_U,{\varphi }_C,{\varphi }_D,x_U,x_C,x_f,x_D\right)$ 。

Figure 1. Three neighboring mesh points and the mesh face

图1. 三个相邻的节点及单元边界

根据Leonard [13] 中的思想，由正则化变量定义：

$\tilde{\varphi }=\frac{\varphi -{\varphi }_U}{{\varphi }_D-{\varphi }_U}$ , $\tilde{x}=\frac{x-x_U}{x_D-x_U}$

正则化前后的变量关系分别如图2，图3所示：

Figure 2. The relationship between variables before the regularization

图2. 正则化前的变量关系

Figure 3. The relationship between variables after the regularization

图3. 正则化后的变量关系

通过正则化公式可以得到 ${\tilde{\varphi }}_U=0,{\tilde{\varphi }}_D=1,{\tilde{x}}_U=0,{\tilde{x}}_D=1$ 。从而可以得到 ${\tilde{\varphi }}_f=f\left({\tilde{\varphi }}_C,{\tilde{x}}_C,{\tilde{x}}_f\right)$ ，我们可以看到 ${\tilde{\varphi }}_f$ 的值是依赖于 ${\tilde{\varphi }}_C,{\tilde{x}}_C$ 和 ${\tilde{x}}_f$ 的值。

下面是一些格式在非结构网格下的线性对流格式。如下表1所示：

Table 1. The linear convection schemes on triangular meshes

表1. 非结构网格下的线性对流格式

2.2. 对流项有界性准则

在数值求解的过程中，在对于对流项的处理过程，既要求具有较高的精度，又要求对于对流项的离散满足对流有界性。因此M. S. Darwish和F. H. Moukalled提出了NVSF方法，并提出了非均匀网格下的高分辨率有限体积格式，并指出由Gaskell和Lau [14] 提出的对流有界性准则CBC对于非结构网格是适用的。在非一致网格下，CBC准则可以写为：

(1)

关于精度问题，Leonard提出在结构网格上利用NVF得出的格式是满足二阶精度的。一个格式为二阶精度的充分必要条件是其表达式过点 $Q\left(0.5,0.75\right)$ ，而对于非一致网格NVSF的情况，所构造的格式是二阶精度的充分必要条件是其表达式过点 $Q\left({\tilde{x}}_C,{\tilde{x}}_f\right)$ 。图4为几种格式在结构网格下与非一致网格下的正则变量图。

Figure 4. Normalized Variable Diagram (NVD) for several linear schemes formulated using NVSF

图4. 利用NVSF构造的几种线性格式的正则变量图(NVD)

2.3. 新格式的构造

我们构造一个新的高阶格式，使该格式满足CBC准则，我们已经知道，在一致网格下，要想使新格式具有二阶精度必须经过点 $\left(0.5,0.75\right)$ 。而在非一致网格下，要想使构造的格式具有二阶精度，必须通过点 $\left({\tilde{x}}_C,{\tilde{x}}_f\right)$ ，我们令

${\tilde{\varphi }}_f=a{\tilde{\varphi }}_C^4+b{\tilde{\varphi }}_C^3+c{\tilde{\varphi }}_C^2+d{\tilde{\varphi }}_C+e$ (2)

此格式需满足：

$\{\begin{array}{l}{\tilde{\varphi }}_f\left(0\right)=0\\ {\tilde{\varphi }}_f\left(1\right)=1\\ {\tilde{\varphi }}_f\left({\tilde{x}}_C\right)={\tilde{x}}_f\\ {{\tilde{\varphi }}^{\prime }}_f\left(1\right)=0\end{array}$ (3)

且我们要构造的新格式与非一致网格下的Quick [15] 格式在 $\left({\tilde{x}}_C,{\tilde{x}}_f\right)$ 处的导数相同，即

${{\tilde{\varphi }}^{\prime }}_f\left({\tilde{x}}_C\right)=\frac{{\tilde{x}}_f\left({\tilde{x}}_f-1\right)}{{\tilde{x}}_C\left({\tilde{x}}_C-1\right)}$ (4)

将(3) (4)带入(2)式中解得

$\{\begin{array}{l}a=\frac{-\left({\tilde{x}}_C^3-3{\tilde{x}}_C^2+3{\tilde{x}}_C{\tilde{x}}_f-{\tilde{x}}_f^2\right)}{{\tilde{x}}_C^2{\left({\tilde{x}}_C-1\right)}^3}\\ b=\frac{-\left(-2{\tilde{x}}_C^4+4{\tilde{x}}_C^3-4{\tilde{x}}_C^2{\tilde{x}}_f+4{\tilde{x}}_C^2+{\tilde{x}}_C{\tilde{x}}_f^2-5{\tilde{x}}_C{\tilde{x}}_f+2{\tilde{x}}_f^2\right)}{{\tilde{x}}_C^2{\left({\tilde{x}}_C-1\right)}^3}\\ c=\frac{-\left({\tilde{x}}_C^5+{\tilde{x}}_C^4-8{\tilde{x}}_C^3+8{\tilde{x}}_C^2{\tilde{x}}_f-2{\tilde{x}}_C{\tilde{x}}_f^2+{\tilde{x}}_C{\tilde{x}}_f-{\tilde{x}}_f^2\right)}{{\tilde{x}}_C^2{\left({\tilde{x}}_C-1\right)}^3}\\ d=\frac{-\left(-2{\tilde{x}}_C^4+4{\tilde{x}}_C^3-4{\tilde{x}}_C{\tilde{x}}_f+{\tilde{x}}_f^2+{\tilde{x}}_f\right)}{{\tilde{x}}_C^2{\left({\tilde{x}}_C-1\right)}^3}\\ e=0\end{array}$

将上式带入(2)式，

可以得到一个非结构网格下的新格式：

当 $0<{\tilde{\varphi }}_C<1$ 时，

$\begin{array}{c}{\tilde{\varphi }}_f=\frac{-\left({\tilde{x}}_C^3-3{\tilde{x}}_C^2+3{\tilde{x}}_C{\tilde{x}}_f-{\tilde{x}}_f^2\right)}{{\tilde{x}}_C^2{\left({\tilde{x}}_C-1\right)}^3}{\tilde{\varphi }}_C^4\\ +\frac{-\left(-2{\tilde{x}}_C^4+4{\tilde{x}}_C^3-4{\tilde{x}}_C^2{\tilde{x}}_f+4{\tilde{x}}_C^2+{\tilde{x}}_C{\tilde{x}}_f^2-5{\tilde{x}}_C{\tilde{x}}_f+2{\tilde{x}}_f^2\right)}{{\tilde{x}}_C^2{\left({\tilde{x}}_C-1\right)}^3}{\tilde{\varphi }}_C^3\\ +\frac{-\left({\tilde{x}}_C^5+{\tilde{x}}_C^4-8{\tilde{x}}_C^3+8{\tilde{x}}_C^2{\tilde{x}}_f-2{\tilde{x}}_C{\tilde{x}}_f^2+{\tilde{x}}_C{\tilde{x}}_f-{\tilde{x}}_f^2\right)}{{\tilde{x}}_C^2{\left({\tilde{x}}_C-1\right)}^3}{\tilde{\varphi }}_C^2\\ +\frac{-\left(-2{\tilde{x}}_C^4+4{\tilde{x}}_C^3-4{\tilde{x}}_C{\tilde{x}}_f+{\tilde{x}}_f^2+{\tilde{x}}_f\right)}{{\tilde{x}}_C^2{\left({\tilde{x}}_C-1\right)}^3}{\tilde{\varphi }}_C\end{array}$

当 ${\tilde{\varphi }}_C\le 0$ 或者 ${\tilde{\varphi }}_C\ge 1$ 时， ${\tilde{\varphi }}_f={\tilde{\varphi }}_C$

所以，非一致网格在退化为结构网格时，点 $\left({\tilde{x}}_C,{\tilde{x}}_f\right)$ 退化为点 $\left(0.5,0.75\right)$ ，则得到结构网格下的新格式为：

$\{\begin{array}{l}{\tilde{\varphi }}_f=-2{\tilde{\varphi }}_C^4+5{\tilde{\varphi }}_C^3-5{\tilde{\varphi }}_C^2+3{\tilde{\varphi }}_C,0\le {\tilde{\varphi }}_C\le 1\\ {\tilde{\varphi }}_f=0,\text{elsewhere}\end{array}$ (5)

3. 二维浅水方程

对于二维浅水方程，我们利用非结构网格下的新格式进行求解。首先考虑二维浅水方程的守恒形式为：

$\frac{\partial U}{\partial t}+\nabla F=0$ (6)

U是守恒变量向量， $F=\left(E,G\right)$ ，E是x方向通量，G是y方向通量。其中

$U=\left(\begin{array}{c}h\\ hu\\ hv\end{array}\right)$ , $E=\left(\begin{array}{c}hu\\ h{u}^2+\frac{g{h}^2}{2}\\ huv\end{array}\right)$ , $G=\left(\begin{array}{c}hv\\ huv\\ h{v}^2+\frac{g{h}^2}{2}\end{array}\right)$

其中，h是水深，u是x方向的平均速度分量，v是y方向的平均速度分量，g是重力加速度。

3.1. 有限体积法数值离散

将计算区域划分为非结构网格，每个三角形单元为控制体 $V_i$ ，把变量定义在三角网格的中心，

在 $V_i$ 上对(6)式积分可得：

${\displaystyle {\int }_{V_i}\left(\frac{\partial U}{\partial t}+\nabla F\right)\text{d}V}=0$

所以可以得到 ${\displaystyle {\int }_{V_i}\frac{\partial U}{\partial t}\text{d}V}+{\displaystyle {\int }_{V_i}\nabla F\text{d}V}=0$

∴ ${\displaystyle {\iint }_A\frac{\partial U}{\partial t}\text{d}A}+{\displaystyle {\iint }_A\nabla F\text{d}A}=0$

∴ ${\displaystyle {\iint }_A\left(\frac{\partial U}{\partial t}+\nabla F\right)\text{d}x\text{d}y}=0$

∴ $\frac{\partial }{\partial t}{\displaystyle {\iint }_{A}U\text{d}A}+{\displaystyle {\iint }_A\nabla F\text{d}A}=0$ (7)

设 $U_i$ 为第i个单元控制体U上的平均值。

$U_i=\frac{1}{A_i}{\displaystyle {\iint }_{A_i}U\text{d}{A}_i},\therefore A_i{U}_i={\displaystyle {\iint }_{A_i}U\text{d}A}$

将其带入(7)式可以得到 $\frac{\partial U_i}{\partial t}=-\frac{1}{A_i}{\displaystyle {\iint }_A\nabla F\text{d}A}$

由Green公式可得：

$\frac{\partial U_i}{\partial t}=-\frac{1}{A_i}{\displaystyle {\oint }_{\Gamma }Fn\text{d}\Gamma }$ (8)

$n$ 为 $\Gamma$ 的单位外法向向量， $\Gamma$ 为控制单元体A的边界。

对(8)式利用Gauss求积公式可得

$\frac{\partial U_i}{\partial t}=-\frac{1}{A_i}{\displaystyle \underset{j=1}{\overset{3}{\sum }}{F}_{ij}{n}_{ij}{l}_{ij}}=-\frac{1}{A_i}{\displaystyle \underset{j=1}{\overset{3}{\sum }}{F}_{ij}\Delta l_{ij}}$

其中， $F_{ij}$ 是通过第i个单元控制体的第j条边上的数值流通量；

$n_{ij}$ 是通过第i个单元控制体的第j条边上的单位外法向量；

$l_{ij}$ 是第i个控制体的第j条边；

$\Delta l_{ij}$ 是第i个单元控制体的第j条边的长度。

通过正则化公式这样就将非一致网格下构造的新格式应用到了三角网格下。

本文采用修正的Roe格式求解界面的数值通量 $F_{ij}$ 。应用经典Roe格式，任意单元体界面通量F可以表

示为 $F_{ij}=h\left({\varphi }_{ij}^L,{\varphi }_{ij}^R\right)$ ，利用Roe [16] 计算数值流通量的式子如下：

$F_{ij}^{Roe}=\frac{1}{2}\left[F\left({\varphi }_{ij}^L\right)+F\left({\varphi }_{ij}^R\right)-\left|a_{ij}\right|\left({\varphi }_{ij}^R-{\varphi }_{ij}^L\right)\right]$ , (9)

其中：

$a_{ij}=\{\begin{array}{l}\frac{F\left({\varphi }_{ij}^R\right)-F\left({\varphi }_{ij}^L\right)}{{\varphi }_{ij}^R-{\varphi }_{ij}^L},{\varphi }^L\ne {\varphi }^R\\ F^{\prime }\left({\varphi }_{ij}^L\right),{\varphi }^L={\varphi }^R\end{array}$ (10)

4. 时间的离散

完成了关于空间导数的离散后，得到了关于时间t的一个常微分方程，即 $\frac{\text{d}\bar{\varphi }}{\text{d}t}=L\left(\varphi \right)$ 。为了抑制非物

理震荡，采用的是三阶Runge-Kutta [17] 格式：

$\{\begin{array}{l}{\varphi }^{\left(1\right)}={\varphi }^n+L\left({\varphi }^{\left(1\right)}\right)\\ {\varphi }^{\left(2\right)}=\frac{3}{4}{\varphi }^n+\frac{2}{3}\left({\varphi }^{\left(1\right)}+L\left({\varphi }^{\left(1\right)}\right)\right)\\ {\varphi }^{\left(n+1\right)}=\frac{1}{3}{\varphi }^n+\frac{2}{3}\left({\varphi }^{\left(n\right)}+L\left({\varphi }^{\left(n\right)}\right)\right)\end{array}$ (11)

5. 数值算例

5.1. 一维线性对流方程

一维线性对流方程

$\{\begin{array}{l}\frac{\partial \varphi }{\partial t}+a\frac{\partial \varphi }{\partial x}=0,x\in \left[a,b\right],t>0\\ \varphi \left(x,0\right)={\varphi }_0\left(x\right)\end{array}$ (12)

在a = 1时，给出不同的初值进行计算。

5.1.1. 情形1

给定以下的初值 $\varphi \left(x,0\right)={\varphi }_0\left(x\right)=\sin \left(2\text{π}x\right)$ ，在区间 $\left[0,1\right]$ 上求解方程(12)，CFL = 0.1，计算时间为0.1，将新构造的高分辨率格式、SMART [9] 格式、HOAB格式、MUSCL格式的数值解进行比较，在计算网格数分别取20、40、80、160、320下，计算格式的 $L_1$ 误差、 $L_2$ 误差以及数值精度阶(如表2)所示，给出格式的精度计算阶公式如下：

${\Vert E\Vert }_p={\left(\Delta x{\displaystyle \underset{i=1}{\overset{N}{\sum }}\left|{\bar{\phi }}_i\left(\text{exact}\right)-{\bar{\phi }}_i\left(\text{computed}\right)\right|}\right)}^{\frac{1}{p}},p=1,2$

$\text{order}=\frac{\ln \frac{E_N}{E_{2N}}}{\ln 2}$

Table 2. Errors and orders for several selected schemes

表2. 格式误差与数值精度阶对比

由表2可知，对于光滑解问题，三种格式的 $L_1,L_2$ 误差随着剖分网格的加密在不断的减小而且结构网格下的新格式可以达到理论上二阶精度。

5.1.2. 情形2

对上述方程，选取间断初值

$\varphi \left(x,0\right)=\{\begin{array}{l}1,x\le 0\\ 0,x\ge 0\end{array}$

选取的计算区域为 $\left[-1,-1\right]$ ，将计算区域等分为800份，计算时间为T = 0.2，由图5可知，新构造的格式在计算间断解时有很好的逼近效果，有效地抑制了非物理震荡。

5.1.3. 情形3

当我们选取如下的形如W信号的非光滑初值时

$\varphi \left(x,0\right)=\{\begin{array}{l}1,0\le x\le 0.2\\ 4x-0.6,0.2<x<0.4\\ -4x+2.6,0.4<x\le 0.6\\ 1,0.6<x\le 0.8\\ 0,\text{elsewhere}\end{array}$

Figure 5. Comparison of numerical and exact results for the new condition

图5. 新格式下T = 0.2时的数值解与精确解

在计算求解过程中，选取的计算区域为 $\left[-0.5,2\right]$ ，网格的剖分份数为320，CFL = 0.1，T = 0.1，下图6给出的是新格式在上述条件下的数值解与准确解的对比，从图上我们可以看到，新格式能很好的计算光滑解，在间断处和大梯度处的近似效果比较好，能有效的抑制非物理震荡。

Figure 6. Comparison of numerical and exact results for the non-smooth initial condition

图6. 在非光滑初值下，新格式的数值解与准确解

5.2. 二维线性对流方程

下面利用非结构网格下的新格式对二维浅水方程进行求解，

$\frac{\partial U}{\partial t}+\frac{\partial E}{\partial x}+\frac{\partial G}{\partial y}=0$

其中， $U=\left(\begin{array}{c}h\\ hu\\ hv\end{array}\right)$ ， $E=\left(\begin{array}{c}hu\\ h{u}^2+\frac{g{h}^2}{2}\\ huv\end{array}\right)$ ， $G=\left(\begin{array}{c}hv\\ huv\\ h{v}^2+\frac{g{h}^2}{2}\end{array}\right)$

U是守恒变量向量，E是x方向通量，G是y方向通量。

我们选取求解区域为 $\left[0,50\right]\times \left[0,50\right]$ ，利用Easy Mesh对求解区域进行三角形网格剖分，图7为三角剖分图。在这种剖分下得到14,770个单元三角形，7546个顶点，22,315条边。从而可以得到二维浅水方程的初始状态的2D和3D图像，如图8所示。利用新构造的格式可以得到数值解，如图9所示，显示了t = 0.49时的数值解的2D和3D图像。由图像可以知道，新的格式在逼近准确解时可以达到很好的效果。

Figure 7. Diagram of triangulation of Easy Mesh

图7. Easy Mesh的三角形剖分示意图

(a) (b)

Figure 8. (a) 2D exact solutions of shallow water equations; (b) 3D exact solutions of shallow water equations

图8. (a) 二维浅水方程的初始状态2D图像；(b) 二维浅水方程的初始状态的3D图像

(a) (b)

Figure 9. (a) 2D numerical solutions of shallow water equations; (b) 3D numerical solutions of shallow water equations

图9. (a) 二维浅水方程t = 0.49时的数值解的2D图；(b) 二维浅水方程t = 0.49时的数值解的3D图

6. 结论

本文是基于三角网格和CBC准则，再结合非一致网格下QUICK格式的导数构造了一个新的NVSF高分辨率格式，用此高分辨率格式对对流项进行离散，时间项采用三阶Runge-Kutta格式进行离散。通过算例的真解与数值解的比较，可以看出所构造的有限体积的高分辨率格式适用于非结构网格，并且将此格式退化为结构网格下的格式也具有很好的精度。

基金项目

本文由内蒙古自然科学基金项目(2015MS0101)和内蒙古自治区人才开发基金项目(12000-1300020240)支持。