基于移动网格的HLL格式求解浅水波方程

Solving Shallow Water Wave Equation with HLL Scheme Based on Moving Grid

摘要: 本文提出了一种基于移动网格的HLL格式来求解浅水波方程。移动网格法基于迭代过程，在每次迭代中利用等分布原理对网格进行重新分布，然后利用保守插值公式对其数值解进行更新，该方法的主要思想是保持每个重新分布步骤下的数值解的质量守恒。HLL格式是具有良好鲁棒性的数值通量，可消除红斑现象。时间方向采用三阶强稳定龙格–库塔方法进行推进，通过数值结果的对比发现基于移动网格的HLL格式具有分辨率高的良好特性。

Abstract: In this paper, an HLL scheme based on moving grid is proposed to solve shallow water wave equations. Based on the iterative process, the moving mesh method redistributes the mesh in each iteration by using the equal distribution principle, and then updates its numerical solution by using the conservative interpolation formula, the main idea of this method is to preserve the mass conservation of the numerical solution under each redistribution step. HLL scheme is a numerical flux with good robustness and can eliminate the carbuncle. The time direction is advanced by the third- order strongly stable Runge-Kutta method, and the comparison of numerical results shows that the HLL scheme based on the moving grid has good resolution.

1. 引言

计算流体力学是指对物理现象进行数值计算时所采用的数值方法，如热传导问题和流体流动问题，然后将数值计算结果显示在图形中。在许多问题中，一些物理量或边界位置在很小的区域内会发生剧烈的变化，自适应网格方法在固体和流体动力学、燃烧、传热、材料科学等各种物理和工程领域都有重要应用，这些领域的物理现象在相当局部的区域，如冲击波、边界层、爆轰波等，发展出动态单一或近似单一的解，这些物理问题的数值研究可能需要在物理区域的一小部分上进行极其精细的网格，以解决大的解变化。在多维情况下，开发有效、鲁棒的自适应网格方法成为必要，成功地实现自适应策略可以提高数值逼近的精度，降低计算成本。Harten和Hyman [1] 开始了这方面最早的研究，经过他们的工作，许多其他移动网格方法的双曲型问题已经被提出且取得显著成果 [2] [3] [4]。在1976年，Yanenko首次提出了一类边界区域大变形问题的求解方法。1983年，Harten提出了一维双曲守恒律的移动网格方法，在每个时间步采用自适应的方法提高了激波和接触间断面的分辨率。从那时起，许多关于双曲线问题的文献被提出。1997年，Li提出将移动网格法和迎风格式相结合来求解含时偏微分方程。该方法保留了r法和h法的优点。1998年，Liu提出了一种基于单元体积的自适应网格方法。该方法应用于欧拉流动。实验结果表明，该方法是获得大梯度(激波等)流场特征清晰分辨率的有效方法。2002年，Azarenok将基于调和映射的变分方法应用于气体动力学的双曲线问题，该方法保持了简单的网格结构，并且所需的计算机代价较小。2010年，钱提出了基于二维欧拉方程移动网格径向基核函数的ENO有限体积格式，减少了计算时间，提高了计算精度，保证了格式的有效性。2017年，Cheng [5] 提出了一种基于移动网格的熵稳定格式用来解决双曲守恒律方程问题，证明该算法更有利于求解大多数区域光滑的问题。2020年，郑素佩，王令等 [6] [7] 针对二维浅水波方程，提出了基于移动网格的熵稳定格式，随后又提出了移动网格旋转通量法，有效提高了浅水波方程数值求解格式的分辨率。

对双曲守恒律方程通量函数的离散方向，Levy等研究了多维控制方程的旋转黎曼求解器，为了更好地捕捉激波，同时保留耗散通量的鲁棒性，研究者采用混合格式。Nishikawa和Kitamura [8] 为了消除激波不稳定性将Roe格式和HLL格式结合起来，其中，HLL格式可以抑制激波不稳定性，Roe格式可以避免过度耗散。Christian Klingenberg [9] 针对具有非平坦底部地形的一维Saint-Venant浅水方程组，构造了新的HLL型移动水守恒迎风格式，通过数值算例验证了所设计的一阶和二阶格式的良好平衡性质以及精确捕捉动水稳态小扰动的能力。贾豆，郑素佩等 [10] 将熵稳定格式与HLL格式结合提出了一种旋转混合格式，用于求解二维欧拉方程，此格式具有数值稳定、分辨率高等优点。

本文对基于移动网格 [4] 的HLL格式进行了研究。介绍了二维浅水波方程、Runge-Kutta法，HLL法等基础知识，通过若干数值算例验证了该格式的性能。

2. 控制方程

二维浅水波方程 [11] 为

$U_t+F{\left(U\right)}_x+G{\left(U\right)}_y=-S,$ (1)

其中，

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

水深h，重力加速度g，x方向上的水速u，y方向上的水速v，源项S。其中S可分为摩擦项和倾斜项，

$S=\left[\begin{array}{c}0\\ gh\frac{\partial z_b}{\partial x}+gh{S}_{f_x}\\ gh\frac{\partial z_b}{\partial y}+gh{S}_{f_y}\end{array}\right],$

$z_b$ 表示底部地势函数， $gh\left(\partial z_b/\partial x\right)$ 和 $gh\left(\partial z_b/\partial y\right)$ 表示水下底部作用力沿x方向和y方向上的分力， $gh{S}_{f_x}$ 和 $gh{S}_{f_y}$ 为水下底面摩擦力沿x方向和y方向的量， $S_{f_x}$ 和 $S_{f_y}$ 为x方向和y方向上的摩阻比率，有

$S_{f_x}=\frac{u}{K^2}\frac{{\left(u^2+v^2\right)}^{1/2}}{h^{4/3}},S_{f_y}=\frac{v}{K^2}\frac{{\left(u^2+v^2\right)}^{1/2}}{h^{4/3}},$

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

3. 自适应移动网格法

自适应移动网格法是基于变分原理的，二维双曲问题自适应网格算法的求解过程具体如下：

步骤1：给定物理域 ${\Omega }_p$ 的一个初始分区 ${\stackrel{\to }{z}}_{j,k}^{\left[0\right]}=\left(x_{j,k}^{\left[0\right]},y_{j,k}^{\left[0\right]}\right)=\left(x_{j,k},y_{j,k}\right)$ 和逻辑域 ${\Omega }_c$ 的一个固定分区，通过单元平均初始数据 $u\left(x,y,0\right)$ 来计算控制体积 $A_{j+1/2,k+1/2}$ 上的网格值 $u_{j+1/2,k+1/2}^{\left[0\right]}$。

步骤2：1) 移动网格 ${\stackrel{\to }{z}}_{j,k}^{\left[v\right]}=\left\{\left(x_{j,k}^{\left[v\right]},y_{j,k}^{\left[v\right]}\right)\right\}$ 到 ${\stackrel{\to }{z}}_{j,k}^{\left[v+1\right]}=\left\{\left(x_{j,k}^{\left[v+1\right]},y_{j,k}^{\left[v+1\right]}\right)\right\}$ 通过Gauss-Seidel迭代法，即

$\begin{array}{l}{\alpha }_{j+\frac{1}{2},k}\left({\stackrel{\to }{z}}_{j+1,k}^{\left[v\right]}-{\stackrel{\to }{z}}_{j,k}^{\left[v+1\right]}\right)-{\alpha }_{j-\frac{1}{2},k}\left({\stackrel{\to }{z}}_{j,k}^{\left[v+1\right]}-{\stackrel{\to }{z}}_{j-1,k}^{\left[v+1\right]}\right)\\ +{\beta }_{j,k+1/2}\left({\stackrel{\to }{z}}_{j,k+1}^{\left[v\right]}-{\stackrel{\to }{z}}_{j-1,k}^{\left[v+1\right]}\right)-{\beta }_{j,k+1/2}\left({\stackrel{\to }{z}}_{j,k}^{\left[v+1\right]}-{\stackrel{\to }{z}}_{j,k-1}^{\left[v+1\right]}\right)=0\end{array}$

其中

${\alpha }_{j\pm \frac{1}{2},k}=w\left(u_{j\pm \frac{1}{2},k}^{\left[v\right]}\right)=w\left(\frac{1}{2}\left(u_{j\pm \frac{1}{2},k+\frac{1}{2}}^{\left[v\right]}+u_{j\pm \frac{1}{2},k-\frac{1}{2}}^{\left[v\right]}\right)\right)$

${\beta }_{j,k\pm \frac{1}{2}}=w\left(u_{j,k\pm \frac{1}{2}}^{\left[v\right]}\right)=w\left(\frac{1}{2}\left(u_{j+\frac{1}{2},k\pm \frac{1}{2}}^{\left[v\right]}+u_{j-\frac{1}{2},k\pm \frac{1}{2}}^{\left[v\right]}\right)\right)$

2) 新网格上的数值解 $\left\{u_{j+1/2,k+1/2}^{\left[v+1\right]}\right\}$ 通过守恒插值进行更新；

3) 重复迭代过程1) 和2) 到指定的迭代步数(3~5步)或直到 $\Vert {\stackrel{\to }{z}}^{\left[v+1\right]}-{\stackrel{\to }{z}}^{\left[v\right]}\Vert \le \epsilon$ ；

步骤3：利用有限体积法在新网格 $\left\{\left(x_{j,k}^{\left[v+1\right]},y_{j,k}^{\left[v+1\right]}\right)\right\}$ 上演化方程，得到时间层 $t_{n+1}$ 上的数值近似 $u_{j+1/2,k+1/2}^{\left[n+1\right]}$ ；

步骤4：如果 $t_{n+1}\le T$，令 $u_{j+1/2,k+1/2}^{\left[0\right]}=u_{j+1/2,k+1/2}^{n+1}$， $\left(x_{j,k}^{\left[0\right]},y_{j,k}^{\left[0\right]}\right)=\left(x_{j,k}^{\left[v+1\right]},y_{j,k}^{\left[v+1\right]}\right)$，返回步骤2直到计算结束。

4. HLL通量格式

在精确Riemann求解器的基础上，Harten等人又提出了HLL通量格式 [12]，该格式是将两种波分离成三个常数状态，HLL格式鲁棒性良好，可以消除红斑现象。HLL格式具体表达式为

$H^{HLL}=\frac{S_R^{+}{H}_k\left(U_L\right)-S_L^{-}{H}_k\left(U_R\right)}{S_R^{+}-S_L^{-}}+\frac{S_R^{+}{S}_L^{-}}{S_R^{+}-S_L^{-}}\Delta U,$ (2)

$H_k=\left(F,G\right)\cdot n\left(k=x,y\right),n=\left(n_x,n_y\right),\Delta U=U_R-U_L,$

$S_R^{+}=\max \left(0,S_R\right),S_L^{-}=\min \left(0,S_L\right),$

$S_L,S_R$ 分别表示左波速和右波速，它们的值分别为Jacobi矩阵 ${\hat{A}}_k=\frac{\partial H_k}{\partial U}$ 特征值的最小值和最大值，即

$S_L=\min \left({\left(q_n\right)}_L-c_L,{\hat{q}}_n-\hat{c}\right),S_R=\max \left({\left(q_n\right)}_R+c_R,{\hat{q}}_n+\hat{c}\right),$

${\left(q_n\right)}_{L,R}={\left(u,v\right)}_{L,R}\cdot n,{\hat{q}}_n=\left(\hat{u},\hat{v}\right)\cdot n,{\hat{q}}_t=\left(\hat{u},\hat{v}\right)\cdot \left(-n_x,n_y\right).$

相对于Godunov方法，给出如下HLL通量，

$H^{\text{HLL}}=\{\begin{array}{ll}{F}_L\cdot n_x+G_L\cdot n_y\hfill & 0\le S_L,\hfill \\ \frac{S_R^{+}{H}_k\left(U_L\right)-S_L^{-}{H}_k\left(U_R\right)}{S_R^{+}-S_L^{-}}+\frac{S_R^{+}{S}_L^{-}}{S_R^{+}-S_L^{-}}\Delta U,\hfill & S_L\le 0\le S_R,\hfill \\ F_R\cdot n_x+G_R\cdot n_y,\hfill & 0\ge S_R.\hfill \end{array}$ (3)

时间方向上进行离散，采用三阶强稳定Runge-Kutta法 [13]，有

$\begin{array}{l}{U}^{\left(1\right)}=U^n+\Delta tL\left(U^n\right),\\ U^{\left(2\right)}=\frac{3}{4}{U}^n+\frac{1}{4}{U}^{\left(1\right)}+\frac{\Delta t}{4}L\left(U^{\left(1\right)}\right),\\ U^{n+1}=\frac{1}{3}{U}^n+\frac{2}{3}{U}^{\left(2\right)}+\frac{2}{3}\Delta tL\left(U^{\left(2\right)}\right).\end{array}$ (4)

5. 数值算例

本节采用所构造的格式求解数值算例，取网格数均为 $50\times 50$，并对所得结果进行分析、比较。

例1 二维圆柱溃坝问题

在计算域 $\left[0,50\text{m}\right]\times \left[0,50\text{m}\right]$ 上，方程(1)当 $S=0$ 时，满足如下条件，

$\begin{array}{l}h\left(x,y,0\right)=\{\begin{array}{l}10,\sqrt{{\left(x-25\right)}^2+{\left(y-25\right)}^2}<11\\ 1,\text{others}.\end{array}\\ u\left(x,y,0\right)=v\left(x,y,0\right)=0,z_b\left(x,y\right)=0.\end{array}$

其中 $CFL=0.45$，时间 $t=1$， $g=0.98$，控制参数 $\alpha =0.835$，边界条件采用周期性条件。图1为圆柱溃坝问题随时间变化的网格演化图；图2为密度等值线图，图2左侧图是在固定网格中由HLL旋转通量格式得到的数值结果，图2右侧图是基于移动网格下的HLL通量得到的结果。可以观察到两种格式均可以捕捉到激波，而基于移动网格下的HLL通量格式得到的数值结果有更高的分辨率。

Figure 1. Grid evolution

图1. 网格演化图

(a) (b)

Figure 2. The cylindrical dam-break density contour. (a) The fixed grid diagram; (b) The moving grid diagram

图2. 圆柱溃坝密度等值线图。(a) 固定网格图；(b) 移动网格图

例2 二维圆形溃坝问题

在区域 $\left[-1,1\text{m}\right]\times \left[-1,1\text{m}\right]$ 上，方程(1)当 $S=0$ 时，满足如下条件，

$\begin{array}{l}h\left(x,y,0\right)=\{\begin{array}{l}2,\sqrt{x^2+y^2}<0.5\\ 1,\text{others}.\end{array}\\ u\left(x,y,0\right)=v\left(x,y,0\right)=0,z_b\left(x,y\right)=0.\end{array}$

其中 $CFL=0.25$，时间 $t=0.2$， $g=0.98$，控制参数 $\alpha =0.1$，边界条件采用周期性条件。图3为圆形溃坝问题随时间变化的网格演化图；图4为密度等值线图。图4左侧图是在固定网格中由HLL旋转通量格式得到的数值结果，图4右侧图是基于移动网格下的HLL通量得到的结果，通过对比可看到基于移动网格下的HLL通量格式得到的数值结果过渡带更窄且无振荡，分辨率更高。

Figure 3. Grid evolution

图3. 网格演化图

(a) (b)

Figure 4. The circular dam-break density contour. (a) The fixed grid diagram; (b) The moving grid diagram

图4. 圆形溃坝密度等值线图。(a) 固定网格图；(b) 移动网格图

例3 二维潮汐模拟问题

方程(1)在区域 $\left[-2,2\text{m}\right]\times \left[-2,2\text{m}\right]$ 上，初始条件与底部地势函数分别为

$h\left(x,y,0\right)=z_b\left(x,y,0\right)+\zeta \left(x,y,0\right),$

$u\left(x,y,0\right)=v\left(x,y,0\right)=0$

$z_b\left(x,y\right)=1+0.01\cos \left(\pi x/2\right)\cos \left(\pi y/2\right)$

其中初始水深表达式为： $\zeta \left(x,y,0\right)=0.1\exp \left(-{10}^{-16}{R}^2\right)$， $R=\sqrt{x^2+y^2}$。

取 $CFL=0.45$，时间 $t=0.3$， $g=0.98$，控制参数 $\alpha =0.935$，边界条件采用周期性条件。图5表示的是二维潮汐问题密度等值线图，网格演化效果不明显，通过对比可以发现两种格式没有明显差别，均对称性良好，结构清晰。

(a) (b)

Figure 5. The density contour of tidal. (a) The fixed grid diagram; (b) The moving grid diagram

图5. 潮汐密度等值线图。(a) 固定网格图；(b) 移动网格图

6. 总结

本文构造了一种基于自适应移动网格的HLL通量格式求解二维浅水波方程。得到的通量函数对激波不稳定具有鲁棒性，分辨率较高。数值实验表明，对于带源项的浅水波方程，两者相比没有明显差别，新算法对称性良好，结构清晰；而对于源项 $S=0$ 的方程，新格式分辨率高，鲁棒性好，具有很好的激波捕捉能力。