非静压模型SWASH学习（2）——潮汐波模拟算例（Tidal wave flow over an irregular bed）

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.2节。
本例考虑了一个不规则地形上的潮汐波，并将数值模拟结果与Bermudez和Vazquez(1994)³ 推导的渐近解作了比较。这个算例能够测试模型是否能处理底坡源项，从而精确模拟水波并避免数值误差引起的波动。
算例的计算域为一个一维的区域（x-方向），其左端有波高4m、波周期12h的潮汐波输入。
以上涉及的数据可下载于传送门。若要运行此算例，请将此下载的压缩文件解压，将含有.sws、 .bnd和 .bot后缀的所有文件解压至swash.exe所在的目录下。

网格、地形及入射波配置

计算区域长 L = 14km，且网格间距 Δx = 280m。其静水深为：
$$d(x)=50.5-\frac{40x}{L}-10sin[\pi (\frac{4x}{L}-\frac{1}{2})]$$
通过上述信息，可以制作地形数据文件 .bot。

入射波从 x = 0 处进入计算域，其变化过程满足：
$$\zeta {|}_{x=0}=4-4sin[2\pi (\frac{t}{43200}+\frac{1}{4})]$$
右侧边界设置为：
$$u{|}_{x=L}=0$$
根据Bermudez和Vazquez(1994)³的研究，该问题的渐进解为：
$$\begin{gathered} \zeta (t)=4-4sin[2\pi (\frac{t}{43200}+\frac{1}{4})] \\ u(x,t)=\frac{\pi (x-L)}{5400(\zeta (t)+d(x))}cos[2\pi (\frac{t}{43200}+\frac{1}{4})] \end{gathered}$$
上式也说明了计算域中的水位仅是时间的函数。

参数设置

网格与地形

首先，设定模型计算域是垂向二维的，即水平向仅有一个维度（ONED）。且模拟的水流不是恒定的（DYN）。

MODE DYN ONED

将计算域划分为矩形网格。计算域总长为 14000.0m，一共划分为50个网格，即网格间距 14000/50 = 280m。

CGRID 0. 0. 0. 14000. 0. 50 0

读取地形数据。地形数据的分辨率设定为10.0m，一共包含1400个网格，正好覆盖了计算域。之后从 'a81tidal.bot’文件中读入数据。

INPGRID BOTTOM 0. 0. 0. 1400 0 10. 0.
READINP BOTTOM 1. 'a81tidal.bot' 1 0 FREE

注：SWASH模型默认读入的BOTTOM数据定义在网格边界点上，而非网格中心点；所以， 'a81tidal.bot’所包含的数字个数为 1400+1 = 1401 个。
由于本算例不考虑垂向

初始条件与边界条件

在初始时刻，模型中的水位值设定为0，所有速度值也设定为0。

INIT zero

之后设定左右边界的边界条件。 首先是设定左边界（W）的条件，即读入入射波的水位变化（BTYPE WLEV）时间序列（SERIES）。这些数据存在于 ‘a81tid01.bnd’ 文件中。注意 最后的参数 itmopt = 7 表示一种时间序列的输入格式。对于itmopt = 7，时间序列文件 ‘a81tid01.bnd’ 内数据需要分两列，第一列为时间数据，写作形如 000107.000 的格式，其等同于 00：01：07.000；第二列为对应的水位数据，单位为m。
右侧边界的速度将保持为0，即右边界的流量（DISCHarge）\通量将保持固定值（CONstant）0.0。

BOU SIDE W CCW BTYPE WLEV  CON SERIES 'a81tid01.bnd' 7
BOU SIDE E CCW BTYPE DISCH CON 0.

物理参数

FRIC CONSTANT 0.
VISC 0.

在本算例中，底面边界通过FRIC指令设定。上述命令表示底面的Manning（默认）系数 cf = 0.0，且cf为固定值（CONSTANT）；即底面为无摩擦的边界。

同时设定水体粘度（VISC）为0.0，即本算例考虑无粘性流动。

数值方法

参数文件中给定了对流项的处理方法，即不采用（NONE）迎风格式（UPW ），而采用标准的中心差分格式。

DISCRET UPW  NONE

此外，对于定义速度的点（速度定义在网格边界上），其水深采用一阶线性（FIRST）插值得到（一般来说，水深、水位值定义在网格中心处）。

DISCRET CORR FIRST

输出控制

本模拟通过如下方式输出结果，首先定义计算域中的一条线，然后在这条线上生成等距的点，并将这些点处的水流数据输出。
QUANTity指令先定义了输出结果的点的位置为一条线上的等距点（DIST），且输出数据的小数点后最多保留位数 HEXP = 10。接着，GROUP指令定义了一个名为LINE的数据包，它包括了第1个到第50个网格的数据。
最后，TABLE指令将数据以一个表格的方式输出，这个表格有文件头（HEAD），即该文件的第一行并非是数据。文件输出了从模拟时间为2：05：52时刻开始的每隔1分钟时刻，各输出点到初始位置的距离（DIST）、底面的静水深（BOTL）、水位（WATL）和速度（VEL）值。

QUANT DIST HEXP 10.
GROUP 'LINE' 1 50  1 1
TABLE 'LINE' HEAD 'a81tid01.tbl' DIST BOTL WATL VEL OUTPUT 020552.000 1 MIN

最后输出文件的头部如下所示：

%
%
% Run:A81   Table:LINE              SWASH version: 7.01
%
%       Dist          Botlev        Watlev      X-vel         Y-vel      
%       [m]           [m]           [m]           [m/s]         [m/s]    
%
      0.00000       60.5000        2.1814      0.114160      0.000000 
    280.00000       59.3858        2.1809      0.114440      0.000000 
    560.00000       57.6631        2.1803      0.115867      0.000000 
    840.00000       55.3897        2.1798      0.118448      0.000000 
   1120.00000       52.6583        2.1791      0.122176      0.000000 
   1400.00000       49.5902        2.1785      0.127003      0.000000 
   1680.00000       46.3279        2.1777      0.132814      0.000000 
   1960.00000       43.0262        2.1769      0.139379      0.000000 
   2240.00000       39.8422        2.1761      0.146305      0.000000 
   2520.00000       36.9258        2.1752      0.152997      0.000000 
   2800.00000       34.4098        2.1743      0.158675      0.000000 

此处的DIST等同于某点到左边界的水平距离（单位：m）。

模拟时间

模拟时间是 00：00：00.000 至 2：06：40.000，且时间步长为1秒。

COMPUTE 000000.000 1 SEC 020640.000

模拟结果

将参数文件 a81tid01.sws，以及地形数据文件 a81tidal.bot 和边界时间序列文件 a81tid01.bnd 复制到swash.exe的同一目录下。并在这个目录下，用如下指令运行：

swashrun a81tid01

全部运行结束后我们会得到对应的.tbl文件和.prt文件。其中，.tbl文件即输出的结果数据；而.prt文件是模型运行过程中的输出。
通过脚本mkplot.m将结果数据绘制成点线图，同时该脚本也将存储在exact.u中的渐进解结果一同绘制在点线图中。matlab脚本文件也在压缩包中。
mkplot.m实现了数据的可视化，于是得到 tidal_1.png 文件（图片如下）。

注：途中的红色圆圈表示渐进解的数据，黑色粗线为SWASH模型模拟结果。

---

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.1016/0045-7930(94)90004-3 ↩︎ ↩︎