非静压模型NHWAVE学习（4）——数值方法

NHWAVE模型采用了Godunov型有限体积法与有限差分法相结合的方式进行数值离散。网格划分后，为确定各控制体积间界面处的通量，NHWAVE模型通过构造界面处的黎曼问题，以其精确解或近似解来确定通量。一般来说Godunov型方法具有良好的激波捕捉能力，可以有效地模拟河口海岸工程中常见的波浪、重力流等模拟问题。
本blog主要基于NHWAVE手册中第三章（3. Numerical Method）的内容，简要介绍了NHWAVE模型中各个部分所采用的数值方法。

网格划分及物理变量的定义

对于黎曼问题，当压力和速度被定义在网格的同一位置时，“棋盘现象”就会发生（详见Patankar或是陶文铨的《数值传热学》教材）。因此，大多数现有模型都使用交错网格，即将压力定义在计算网格的中心、速度定义在网格间的界面中心；简单来说，就是将压力点与速度点分开布置。不过，这种传统的交错网格不像同位网格那样，易于应用Godunov型格式。同时，在自由水面应用压力边界条件时，可能会出现无法准确描述顶层网格中心压力值的困难¹。
考虑到上述情况，NHWAVE中采用了一种名为Keller-Box²的离散物理量定义方式（如下图所示）。在这种模式下，速度值被定义在网格的中心（对流扩散方程中的浓度项C也被定义在网格中心），而压力定义在网格垂向界面的中心。
（注：在模型计算中，水深/水位也定义在各个水平网格的中心）

这样一来，位于自由水面出的压力可以得到精确地表示（p = 0 at the free surface）。
还有一点需要说明的是，目前的NWAVE模型之能支持结构化矩形网格，而且在同一方向上，网格间距是一个固定值。

时间积分格式（Runge-Kutta Methods）

在动量方程中，有一个速度关于时间的变化率项，即下图式子的第一项。

对于该项的处理，NHWAVE模型中设置了三个可选项，分别是一阶向前差分、二阶Runge-Kutta法、和三阶Runge-Kutta法。

1. 一阶差分法。即在动量方程的第一项写成如下形式：
   $$\frac{\partial U}{\partial t}=\frac{U^{n+1}-U^n}{\Delta t}$$
   其中，上标n表示当前时间步，n+1表示后一个时间步；并在后面计算F、G、H和源项时，都采用n时间步的物理量。

2. 二阶Runge-Kutta法。将对时间变化率的求解分为以下两个阶段：
   
   和
   
   其中U^*为中间量，而商标（1）和（2）表示计算的两个阶段。在R-K法的每个阶段，模型分别进行了两次计算，一次是考虑除了S_p以外的所有项，第二次是仅考虑了S_p的效应，这相当于先不考虑动水压力计算一次，而后再利用动水压力进行修正；最后，再根据上述求得的两个阶段值，用以下的方式给出加权平均值，以作为下一个时间步的U。
   

3. 三阶Runge-Kutta法。方法与二阶Runge-Kutta法类似，同样是在每一个阶段先不考虑S_ρ计算一次，之后再用S_ρ修正一次；不同的是，在第二阶段计算结束后，模型通过如下方式得到第三阶段的U：
   $$U^{(3)}=\frac{3}{4}{U}^{(1)}+\frac{1}{4}{U}^{(2)}$$
   最后，通过下式得到n+1时间步的U：
   $$U^{n+1}=\frac{1}{3}{U}^{(2)}+\frac{2}{3}{U}^{(3)}$$

一般来说，在同等的网格下，方法的阶数越高，得到数值解的精度也就越高，但消耗的计算资源也就越大。而且，求解S_p时需要迭代求解压力泊松方程，这个步骤正是模型计算中最为耗时的。有研究表明，在NHWAVE模型算例中，求解泊松方程的时间占总耗时的70%以上³。
无论是几阶时间积分格式，在每一阶段，模型都是通过显式求解式得到水位、速度等物理量。在此过程中，模型的时间步长是自适应变化的；时间步长的大小遵循下面的Courant-Friedrichs-Lewy (CFL) 准则：

其中，C表示Courant数。保持格式稳定的必要条件是C<1.0；在一般的算例中，我们常取C=0.5。

空间离散格式

基于TVD方法的空间物理量重构

对于动量方程中的各个通量项（如下图所示），它们都定义在网格的边界上。然而，下列表达式中的速度、水位等各个值并非位于网格的边界上，故我们需要通过一定的方法，以得到边界上各个物理量的值。

NHWAVE模型采用了一种分段线性重构的方法⁴，以得到各个边界处的物理量值。在此，我们假定一个定义在网格 i 中心的物理量U_i；于是在网格内，该物理量可表达为：
$$U=U_i+(x-x_i)\Delta U_i,x\in [x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$$
其中，x表示某一方向的坐标值，ΔU_i表示该网格中物理量的斜率（slope）。上述方法满足TVD条件，经证明可以在解的平滑处能得到二阶精读，且能够高效地捕捉到激波现象。对于ΔU_i，可通过下式确定：
$$\Delta U_i=limter(\frac{U_{i+1}-U_i}{x_{i+1}-x_i},\frac{U_i-U_{i-1}}{x_i-x_{i-1}})$$
其中，limter() 表示限制器函数（slope limiter）。在NHWAVE的源代码中，配置了三种limiter供大家选择，它们分别是
$$\begin{gathered} Minmod:limter(a,b)=max[0,min(a,b)]; \\ vanLeer:limiter(a,b)=\frac{a|b|+|a|b}{|a|+|b|}; \\ superbee:limiter(a,b)=sgn(b)\cdot max[0,min(2|b|,sgn(b)\cdot a),min(|b|,sgn(b)\cdot 2a)] \end{gathered}$$

界面通量的求解

通过上面的方式，我们能得到各个界面上物理量的值；不过，每个界面上得到的值将存在两个状态。
以界面x_i+1/2为例，它的左侧是以x_i为中心的网格，右侧是以x_i+1为中心的网格；在这两个网格中分别使用分段线性重构，我们可以得到：
$$\begin{gathered} U_{i+\frac{1}{2}}^L=U_i+\frac{1}{2}\Delta x_i\Delta U_i \\ {U}_{i+\frac{1}{2}}^R=U_{i+1}-\frac{1}{2}\Delta x_{i+1}\Delta U_{i+1} \end{gathered}$$
其中，商标L和R分别表示该界面的左、右侧状态值。基于该两个状态值，NHWAVE采用近似Riemann问题求解器得到界面通量F；在NHWAVE模型中，配置的求解器有两种，一个是HLL，一个是HLLC，它们的表达式如下：

1. HLL solver
   
   
2. HLLC solver
   
   式中，D表示水深，水深和速度的下表L、R分别表示该网格界面左、右两侧对应网格中心的水深或速度值。
   对于它们之间的差异，可参考文献5⁵。

压力泊松方程

在介绍时间积分时候我们提到，我们需要通过动水压力p来修正不考虑S_p时得出的流场状态，即：

其中，商标 () 表示时间积分Runge-Kutta法的计算阶段。
而为了得到动水压力p，我们需要求解压力泊松方程。在该方程中，只有同一时间步下的物理量和空间变化率项，不存在时间变化率项。为了得到该方程，我们需要将上式带入σ坐标系下的连续性方程

可以化简得到：

利用二阶空间中心有限差分格式对上述方程进行离散化。速度(u, v, w)在网格垂直界面的值通过相邻网格插值得到。所得到的线性方程组有如下形式：


所构成的线性方程组含有不对称的系数矩阵，共有15条对角线；利用高性能HYPRE库上述线性系统进行求解，以得到动压力p。

其它

边界条件

1. 自由水面的水平速度条件（不考虑风应力）
   

2. 自由水面的垂向速度条件；设定虚网格，给出虚拟网格垂直速度，以满足
   

3. 自由水面的压力条件
   

4. 底部的垂向速度条件
   

5. 底部的水平速度条件
   free-slip条件：
   
   考虑摩阻力的边界条件：
   
   其中，u_b表示水平速度的矢量 (u, v)，c_d表示阻力系数.若底面为非滑移的fully rough，c_d 可通过壁面函数（the law of the wall）得到；若给定了底面粗糙高度k_s，则：
   

6. 底部动压力边界条件（Neumann条件）
   
   注：在NHWAVE的landslide模块中，若考虑属下滑坡体，则动压边界条件可线性化为
   

在封闭边界或垂直壁面处，模型一般是应用自由滑移边界条件，即使法向速度和切向应力均为零、法向压力梯度为零。
在考虑波浪输入或入流时，则根据解析解计算结果指定边界处水深和速度。
在横向上，我们还可以应用周期性边界条件。（注：使用并行计算时，与周期性边界垂直的方向上的划分网格数必须是2的幂）
为了方便并行计算，我们在每个边界上使用了两个虚网格。边界条件一般在这两个虚网格中指定。

对干湿点的处理

每个网格的干湿情况都是通过计算出的水深D来判断。首先指定一个小水深D_min，如果D>D_min，则为湿网格（参与计算），反之则为干网格（暂时不参与计算）。
在源代码中，有一变量Mask指示了网格的干湿情况；Mask = 1 表示湿网格，Mask = 0 表示干网格。设定 η_i,j = D_min - h_i,j, h表示静水深。之后通过下列条件判断各个干网格是否要重新激活，变成湿网格：

在干网格中，网格界面处的通量设为零；此外，干网格对应的通量求解步骤中，需要采用以下的修正式：

对于对流-扩散方程的处理

对于上述的对流扩散方程，可采用与动量方程类似的方法离散。不过在求解对流扩散方程是，式中的水深、速度等水动力物理量均为已知值。
在σ坐标系下，对流扩散方程一般有如下形式

注意，其中的ω_s表示沉降速度，当考虑泥沙输运时，该值需要被计算，而没有考虑泥沙时，该值被设定为0；σ_h,σ_v分别表示水平和垂直方向上的施密特数。

参考文献

---

1. Yuan H , Wu C H . An implicit three‐dimensional fully non‐hydrostatic model for free‐surface flows[J]. International Journal for Numerical Methods in Fluids, 2004, 46(7):709-733. ↩︎

2. Stelling G , Zijlema M . An accurate and efficient finite-difference algorithm for non-hydrostatic free-surface flow with application to wave propagation[J]. International Journal for Numerical Methods in Fluids, 2003, 43(1):1-23. ↩︎

3. 时健,郑金海,严以新,等. 河口海岸水动力非静压数学模型研究述评[J]. 河海大学学报（自然科学版）,2017,45(2):167-174. ↩︎

4. Zhou J G , Causon D M , Mingham C G , et al. The Surface Gradient Method for the Treatment of Source Terms in the Shallow-Water Equations[J]. Journal of Computational Physics, 2001, 168( 1):1-25. ↩︎

5. Ma G , Kirby J T , F Shi. Numerical simulation of tsunami waves generated by deformable submarine landslides[J]. Ocean Modelling, 2013, 69(Complete):146-165. ↩︎