非静压模型SWASH学习（3）——水跃模拟算例（Hydraulic jump in an open channel）

SWASH是由Delft大学开发，用于模拟非静压条件下的水动力/波浪运动的数值模型。
与模型原理相关的内容详见以下论文：

1. SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters (Marcel Zijlema, Guus Stelling, Pieter Smit)¹
2. Computation of free surface waves in coastal waters with SWASH on unstructured grids (Marcel Zijlema)²

其中，第2篇论文是描述了SWASH模型的最新功能，即可支持非结构化三角形网格。不过，本blog提及的算法、设置以swash结构化网格的版本为对象，不涉及非机构化网格。
对于SWASH模型，我们从论文1提及的各个验证算例开始，介绍模型配置和参数的含义。
注：有些指令、参数已经在前面的博客中进行讲解了，故之后不会再详细说明；相关内容详见之前的博文。也希望大家能留言，来相互交流！

SWASH主页：https://swash.sourceforge.io/
模型手册：https://swash.sourceforge.io/online_doc/swashuse/swashuse.html

算例简介

本blog所涉及算例详见论文《SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters (Marcel Zijlema, Guus Stelling, Pieter Smit)》¹ 的5.5节。
本算例模拟了水流从急流向缓流转变过程中产生的水跃现象。当水跃发生时，水面会急剧上升，发生湍流混合、能量耗散等。通过此算例，我们评估了SWASH模型预测明渠水流变化及模拟水跃现象的能力。该算例的模型配置同 Zhou and Stansby (1999)³ 文中的数值实验。
模型输入文件及相关数据可下载于【传送门】。若要运行此算例，请将此下载的压缩文件解压，将含有.sws、 .wlv和 .bot后缀的所有文件解压至swash.exe所在的目录下。

模型配置

模型计算域长 L = 30.5m。地形由坡度不同的两部分组成：第一部分是水平的，长14.5 m；第二部分长16.0 m，坡度为0.03（见下图，图中粗实线表示底面）。

水流从计算域的左侧流入，水流的流速 u 满足：
$$Fr=\frac{u}{\sqrt{gh}}=4.65$$
其中，h表示入流水深。

参数设置

网格与地形

MODE DYN ONED
CGRID 0. 0. 0. 30.5 0. 60 0

本算例计算域为一维（ONEDimension），且模拟的水流为非定常流动（DYNamic）。
计算域被均匀划分为矩形网格（CGRID），区域总长30.5m，网格总数为60。

INPGRID BOTTOM 0. 0. 0. 60 0 0.5 0.1
READINP BOTTOM 1. 'a46jump.bot' 1 1 FREE

网格数据在a46jump.bot文件中，文件中共包含61个数据。

初始条件

模型的初始条件设定过程分为两部分。第一部分是导入初始时刻的水面线数据（WLEV），该数据在文件 'a46jump.wlv’中；第二部分是设定初始速度，在本算例中，初始速度将通过Chezy公式从水位推导出来。上述设定方式可以缩短SWASH模型的启动时间，这在准稳态流动条件下（例如模拟河流中的水流）非常实用。

INPGRID WLEV 0. 0. 0. 1 0 30.5 1.
READINP WLEV 1. 'a46jump.wlv' 1 0 FREE
$
INIT steady

首先定义了读入的初始水位场是一个长 30.5m，且被离散为只有一个1个网格的数据。所有的水位数据将被存在网格的边界点上（而非网格中心，所以所需要的数据点总共有两个）。而具体的数据被存储在’a46jump.wlv’文件中。 'a46jump.wlv’文件中的数据如下：

0.06
0.06

这表示初始时的水位为 η = 0.06m。
之后，通过 INITIAL steady 指令启用Chezy公式来推导初始速度场。
在本算例中，除入流边界外，其余位置的初始速度的计算结果为0；而在入流边界处，上述方法将定义一个黎曼不变量，来作为输入值，该值 u_inlet的值通过下式得出：
$$u_{inlet}=u+2\sqrt{gh}=Fr/\sqrt{gh}+2\sqrt{gh}$$

边界条件

BOU SIDE W CCW BTYPE RIEMANN CON 5.105
BOU SIDE E CCW BTYPE OUTFLOW

对于左边界（West），入流的边界通过定义一个黎曼不变量来实现（见上一小节），这个量的值为定值，大小为5.105。
模型的右边界（East）为一个出流边界（OUTFLOW）。在SWASH模型中，出流边界的实现方式为边界处的零水深梯度，即：
$$\frac{\partial h}{\partial x}=0$$

物理参数

FRIC MANNING 0.019
VISC 0.

模型的底摩擦系数（Manning系数）设定为0.019。水体粘度设定为0，即本算例模拟的为无粘性流动。

数值格式

DISCRET UPW  MOM MUSCL
DISCRET CORR FIRST

对于水平项的对流过程，模型采用了动量守恒格式（MOMentum），并在计算通量时采用MUSCL限制器。此外，对于流速计算点的水位和水深，模型采用一阶（FIRST）迎风插值的方法。

注：关于水动力模型中限制器的相关描述和研究，可查看本人的论文 Effects of Different Slope Limiters on Stratified Shear Flow Simulation in a Non-hydrostatic Model (Hu, 2022)⁴。

TIMEI METH IMPL 1.

对于时间离散，模型采用半隐式格式（simi-IMPLicit）。当用了半隐式格式后，水位梯度项和连续方程中采用的速度散度都需要用 θ−method 进行半隐式离散。这种方法在模拟由浮力、潮汐或风驱动的三维环流时特别有效，因为在这些情况下，库朗数常会大于1，但模型仍能提供足够精确的解。
一般来说，需要分别给定θ−method运用在连续方程和水位梯度计算时的θ值，即模型中的thetac和thetas参数。在本算例中，thetac=1.0，即对于连续性方程，模型采用了全隐的Euler格式。（对于水位梯度，thetas默认为0.5，即采用 Crank-Nicolson 格式）。

更多说明详见SWASH手册。

输出控制

QUANT DIST HEXP 10.
GROUP 'LINE' 1 61  1 1
TABLE 'LINE' HEAD 'a46jmp02.tbl' DIST BOTLEV DEPTH WATL VEL

模型将输出所有网格边界点位置上的水平距离（DIST，指到原点的距离）、底高程（BOTLEV）、水深（DEPTH）、水位（WATL）和水平速度（VEL）数据。
由于网格总数为60，故网格边界点有61个，所以GROUP后续的参数为 1 61，这两个数用于表示起始的网格边界序号和末位网格边界序号。
在TABLE指令中，我们没有指定输出数据的模拟时间，即没有设定OUTput参数。此时，模型默认输出最后时刻的流场数据。

计算时间

COMPUTE 000000.000 0.1 SEC 001000.000

模拟的起始时间被设定为 00：00：00.000，时间步长为0.1秒；模拟结束时间为00：10：00.000。

模拟结果

将参数文件 a46jmp02.sws，以及地形数据文件 a46jmp.bot 和初始水位数据文件 a46jump.wlv 复制到swash.exe的同一目录下。并在这个目录下，用如下指令运行：

swashrun a46jmp02

运行完成后，得到结果数据文件 a46jmp02.tbl。
通过脚本mkplot.m将结果数据绘制成点线图。matlab脚本文件也在压缩包中。同时，mkplot.m实现了数据的可视化，于是得到 jump_7.png 文件（图片如下）。

最后，为了验证计算的精确度，我们将模拟结果与Bélanger公式计算结果进行比较。Bélanger公式如下：
$$\frac{h_1}{h_0}=\frac{-1+\sqrt{1+8F{r}^2}}{2}$$
式中，h₀和h₁分别表示水跃的跃前水深与跃后水深；它们的比值取决于上游侧的Froude数。由模拟结果可知，h₀ = 0.1102m（位置x≈5.1m），h₁ = 0.2194m（位置x≈7.625m），则h₁/h₀=1.99；而Bélanger公式的计算结果为（取x≈5.1m处的Froude数1.78）
$$\frac{h_1}{h_0}=\frac{-1+\sqrt{1+8F{r}^2}}{2}=2.07$$
模拟和公式计算值的相对误差约为3.9%，这说明模拟结果的精度是满足一般工程要求的。

---

1. https://doi.org/10.1016/j.coastaleng.2011.05.015 ↩︎ ↩︎

2. https://doi.org/10.1016/j.compfluid.2020.104751 ↩︎

3. https://doi.org/10.1006/jcph.2000.6670 ↩︎

4. https://doi.org/10.3390/jmse10040489 ↩︎