求解二维浅水波方程的旋转混合格式*

郑素佩， 李 霄， 赵青宇， 封建湖

（长安大学 理学院，西安 710064）

引 言

随着计算机技术的不断发展，对于计算流体动力学中更精确的数值计算方法的需求也越来越迫切.浅水波方程是浅水环境下流体运动的一种数学描述，是研究河流、灌溉和海洋等波动问题的一种重要模型，对其的精确求解还存在很大困难.浅水波方程本质上是一类非线性双曲守恒方程，对于双曲守恒律方程数值求解算法的研究一直是一个热点，间断解的存在对该类方程的数值求解格式提出了很高的要求.Tadmor[1]提出了熵稳定格式，为得到稳定的具有物理意义的解提供了一种简单、有效的方法.Roe[2]在自己提出的熵守恒格式的基础上增加了数值黏性项，得到的熵稳定格式能很好地捕捉激波，抑制间断处出现的伪振荡.程晓晗、蔡力等[3]提出了一种基于WENO 重构的熵稳定格式，该格式具有基本无振动性和高分辨率的特点.Liu 等[4]采用SPWENO/WENO-YC3 重构构造了高分辨率熵稳定格式，该格式具有稳定、低耗散的特点.王令、郑素佩[5]针对二维浅水波方程，提出了基于移动网格的熵稳定格式，有效提高了浅水波方程数值求解格式的分辨率.

对双曲守恒律方程通量函数的离散方向，Levy 等[6]研究了多维控制方程的旋转Riemann 求解器，讨论了通量函数沿迎风方向和界面法线方向的离散数值求解算法.Ren 提出了一种旋转Roe 格式[7]，该格式不受激波不稳定性的影响，从理论上分析了旋转通量函数的耗散特征.Zhang 等[8]提出了一种旋转迎风格式求解二维问题，其结果能够明显减少额外的耗散.郑素佩等[9]基于移动网格采用旋转通量求解二维浅水波方程，该算法具有高分辨率的特点.为了更好地捕捉激波，同时保留耗散通量的鲁棒性，研究者们采用了混合格式.Nishikawa 和Kitamura[10]为了消除激波不稳定性将Roe 格式和HLL 格式结合起来，其中，HLL 格式可以抑制激波不稳定性，Roe 格式可以避免过度耗散.刘友琼等[11]基于旋转Riemann 求解器将HLLC 格式与HLL 格式进行特定结合得到了一类混合型数值格式，此格式具有结构简单、数值稳定、分辨率高等优点.贾豆、郑素佩[12]将熵稳定格式与HLL 格式结合提出了一种旋转混合格式，用于求解二维Euler 方程，数值结果分辨率高.

本文对由熵稳定格式和HLL 格式结合得到的旋转混合格式进行了研究，介绍了二维浅水波方程、有限体积法与时间方向离散的Runge-Kutta 法的基础知识，给出了旋转通量混合格式的具体表达式.通过若干数值算例验证了该旋转混合格式的性能，表明该混合格式对于二维浅水波方程数值求解鲁棒性好、分辨率高、具有很好的激波捕捉能力.

1 控制方程

二维浅水波方程[13]为

其中

h表示水深，g表示重力加速度，u,v分 别表示x方 向和y方 向上的水流速度，S表示源项.其中S可分为摩擦项和倾斜项：

zb表示底部地势函数，gh(∂zb/∂x)和gh(∂zb/∂y) 表示水下底部作用力沿x方 向和y方向上的分力，ghSfx和ghSfy为水下底面摩擦力沿x方 向和y方 向的分量，Sfx和Sfy为x方 向和y方向上的摩阻比率，有

K表示摩擦因数，通常情况取K=50.

2 方程离散

本文采用任意四边形单元对计算区域进行离散[9]，即控制体为任意四边形单元.

2.1 旋转不变性[5,14]

对浅水波方程（1）左右两边同时进行积分，得到

其中V表示控制体，A表示V的边界，H=(F,G)表示通量的张量，n为 边界面A的单位外法向量， dA表示一个面积元，H·ndA表示边界A的通量法向分量.

方程（3）的半离散格式如下：

2.2 数值通量

2.2.1 熵稳定格式

二维浅水波方程熵稳定格式可表示为下面的形式[15-19]：

2.2.2 HLL 通量格式

在精确Riemann 求解器的基础上，Harten 等提出了HLL 通量格式[20-23]，该格式是将两种波分离成三个常数状态，具有良好的鲁棒性，能有效地消除红斑现象.HLL 格式的具体表达式为

2.2.3 混合通量格式

2.3 时间离散

采用三阶强稳定Runge-Kutta 法[24]对时间进行离散，有

3 数值算例

本节采用所构造的混合格式对几个浅水波方程数值算例进行求解，取网格数为1 50×150，并对所得结果进行分析、比较.

例1二维圆柱溃坝问题

考虑方程（1）当S=0 时在计算域[ 0,50]×[0,50]上满足如下条件的问题：

其中 C FL数为0.45 ， 时间t=1,g=0.98，采用周期性边界条件.图1～3 分别表示二维圆柱溃坝问题的模拟图、密度等值线图和速度等值线图.图2(a)、3(a)和3(c)是仅由熵稳定旋转通量格式得到的数值结果，图2(b)、3(b)和3(d)是用本文构造的旋转混合通量得到的结果，可以看出两种格式均可以捕捉到激波，而旋转混合格式得到的数值结果具有更高的分辨率.

图1 圆柱溃坝模拟图Fig.1 The simulation diagram for the cylindrical dam

图2 二维圆柱溃坝问题密度等值线图：(a) 非混合格式密度结果；(b) 混合格式密度结果Fig.2 The density contour for the 2D cylindrical dam-break problem：(a) the non-mixed scheme density results; (b) the mixed scheme density results

例2二维区域上的圆形溃坝问题

考虑方程（1）当S=0 时在区域[ −1,1]×[−1,1]上满足如下条件的问题：

其中C FL数为0.25 ， 时间t=0.2,g=0.98，采用周期性边界条件.

图3 二维圆柱溃坝问题速度图：(a) 非混合格式x 方向的速度u；(b) 混合格式x 方向的速度u；(c) 非混合格式y 方向的速度v；(d) 混合格式y 方向的速度vFig.3 The velocity contour for the 2D cylindrical dam-break problem: (a) the non-mixed scheme x direction velocity u results; (b) the mixed scheme x direction velocity u results; (c) the non-mixed scheme y direction velocity v results; (d) the mixed scheme y direction velocity v results

图4～6 分别表示圆形溃坝问题的模拟图、密度等值线图和速度等值线图.图5(a)、6(a)和6(c)是仅由熵稳定旋转通量格式得到的数值结果，图5(b)、6(b)和6(d)是由旋转混合格式得到的结果，通过对比可以看到旋转混合格式得到的数值结果具有更高的分辨率，对称性好，没有出现非物理振荡.

图4 圆形溃坝模拟图Fig.4 The simulation diagram for the circular dam break

图5 二维圆形溃坝问题密度等值线图：(a) 非混合格式密度结果； (b) 混合格式密度结果Fig.5 The density contour for the 2D circular dam-break problem: (a) the non-mixed scheme density results; (b) the mixed scheme density results

图6 二维圆形溃坝问题速度图：(a) 非混合格式x 方向的速度u；(b) 混合格式x 方向的速度u；(c) 非混合格式y 方向的速度v；(d) 混合格式y 方向的速度vFig.6 The velocity contour for the 2D circular dam-break problem: (a) the non-mixed scheme x direction velocity u results; (b) the mixed scheme x direction velocity u results; (c) the non-mixed scheme y direction velocity v results; (d) the mixed scheme y direction velocity v results

例3二维区域上的激波聚焦问题

考虑方程（1）当S=0 时在区域[ −1,1]×[−1,1]上满足如下条件的问题：

其中C FL数为0.6,时 间t=0.2,g=0.98,采用周期性边界条件.图7～9 分别表示二维激波聚焦问题的模拟图、密度等值线图和速度等值线图.图8(a)、9(a)和9(c)是仅由熵稳定旋转通量格式得到的数值结果，图8(b)、9(b)和9(d)是由旋转混合格式得到的结果，可以看出旋转混合格式过渡带更窄且无振荡，具有更高的分辨率.

图7 激波聚焦模拟图Fig.7 The simulation diagram for the shock wave focusing problem

图8 二维激波聚焦问题密度等值线图：(a) 非混合格式密度结果； (b) 混合格式密度结果Fig.8 The density contour for the 2D shock wave focusing problem: (a) the non-mixed scheme density results; (b) the mixed scheme density results

图9 二维激波聚焦问题速度图：(a) 非混合格式x 方向的速度u； (b) 混合格式x 方向的速度u；(c) 非混合格式y 方向的速度v ； (d) 混合格式y 方向的速度vFig.9 The velocity contour for the 2D shock wave focusing problem: (a) the non-mixed scheme x direction velocity u results; (b) the mixed scheme x direction velocity u results; (c) the non-mixed scheme y direction velocity v results; (d) the mixed scheme y direction velocity v results

例4二维潮汐模拟问题

在区域[ −2,2]×[−2,2]上考虑方程(1)，底部地势函数和初始条件为

zb(x,y)=1+0.01cos(πx/2)cos(πy/2),

h(x,y,0)=zb(x,y,0)+ζ(x,y,0),u(x,y,0)=v(x,y,0)=0,

其中初始水深表达式为

取CFL数为0.45， 时间t=0.3,g=0.98, 采用周期性边界条件.图10 为二维潮汐问题的水面高度模拟图，图11为潮汐问题的密度等值线图，图12 为二维潮汐问题的速度图.图11(a)、12(a)和12(c)是仅由熵稳定旋转通量格式得到的数值结果，图11(b)、12(b)和12(d)是由熵稳定和HLL 耦合的旋转混合格式得到的结果，通过对比可以发现两种格式没有明显差别， 均对称性良好，结构清晰.

图10 潮汐问题模拟图Fig.10 The simulation diagram for the tidal problem

图11 二维潮汐问题密度等值线图：(a) 非混合格式密度结果；(b) 混合格式密度结果Fig.11 The density contour map for the 2D tidal problem: (a) the non-mixed scheme density results; (b) the mixed scheme density results

图12 二维潮汐问题速度图：(a) 非混合格式x 方向的速度u；(b) 混合格式x 方向的速度u；(c) 非混合格式y 方向的速度v；(d) 混合格式y 方向的速度vFig.12 The velocity contour for the 2D tidal problem: (a) the non-mixed scheme x direction velocity u results; (b) the mixed scheme x direction velocity v results; (c) the non-mixed scheme y direction velocity v results; (d) the mixed scheme y=direction velocity v results

4 总 结

本文构造了一种求解二维浅水波方程的旋转通量混合格式.利用浅水波方程的旋转不变性将熵稳定格式和HLL 格式耦合得到混合格式进行数值模拟时，对于带源项的浅水波方程，混合旋转通量法与非混合格式相比没有明显差别，对称性良好，结构清晰；对于源项为零的浅水波方程，混合旋转通量法具有更高的分辨率.验证可得，新的旋转混合格式鲁棒性好、分辨率高、具有很好的激波捕捉能力.