非静压模型SWASH学习（1）——线性波模拟算例（Linear progressive waves through a flume）

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

算例简介

本算例模拟了一个水槽中正在行进的线性波，以验证模型精确解析频散的能力。这个算例是验证线性频散精度的一个非常常见的测试算例。
首先，我们建立一个垂向二维（有x-和z-方向）的数值水槽，其左端有波长为 λ 的入射波，水槽总长为 36λ。基于线性波动理论，在水槽的波入口处施加了弱反射双曲余弦波的垂直速度分布。在水槽右端应用了宽度为 36λ 的海绵层和 Sommerfeld 辐射条件，以最大限度地减少波浪反射。
本算例一共考虑了 kd = π、2π、3π、4π 和 5π 五种情况，而且满足 A₀/d = 0.001；其中A₀表示波振幅，d表示静水深。此外，波周期定义为T₀。

参数文件解读

为运行此算例，我们需要一个参数文件（.sws）及一个水深文件（.bot）。参数文件中包含了模型运行所需的所有参数，以及数值方法的选择、网格设置和边界条件设置等。水深文件给定了计算域的地形数据。以上涉及的数据可下载于传送门。若要运行此算例，请将此下载的压缩文件解压，将含有.sws和.bot后缀的所有文件解压至swash.exe所在的目录下。
在本节中，我们来以a14prw01.sws为例，借助模型手册来解读参数文件。

参数文件结构

该参数文件可分为文件头（HEADING）、输入（MODEL INPUT）和输出（OUTPUT REQUESTS）这三部分，最终以STOP指令表示结束。住：在.sws文件中，$ 或 ! 符号表示注释，即以 $ 或 ! 开头的行不会被读取，仅有注释作用。
其中文件头包含了算例的名称；输入部分包含了网格、初始条件、边界条件、数值方法等；输出部分定义了文件格式、输出频率等。

文件头（HEADING）

文件头的格式如下：

PROJect ’name’ ’nr’
’title1’
’title2’
’title3’

其中，name填的是模拟算例的名称，最长16个字符；nr表示运行标识的总数，以区分name相同的各个算例；下面的title1、title2… 都是一个最多含72个字符的字符串，该字符串会出现在程序的输出中文件中。注：PROJect指令在输入时可以省略其中的几个小写字母，而简写成PROJ。
本算例的文件头是：

PROJ 'A14prw01' 'A14'

它表示了这是算例A14prw01中名称为A14的数值模拟。

计算区域和网格

模型的计算网格和区域由以下指令设定：

$
MODE DYN ONED
$
CGRID 0. 0. 0. 240. 0. 360 0
$
VERT 2
$
INPGRID BOTTOM 0. 0. 0. 1 0 240. 1.
READINP BOTTOM 10. 'a14prwav.bot' 1 0 FREE
$

第一条指令 MODE 表示了模拟的类型，该指令有两个参数；

第一个参数可填 NONS（或DYN），表示为非平稳流场的模拟；若不填，则表示稳态流动的模拟。第二个参数选填 TWOD和ONED，分别表示水平方向上的二维和一维模拟。

第二的指令CGRID用于计算网格的生成。在本算例中，我们采用结构化的矩形网格（REGular）。

一共需要七个参数，这七个参数的含义可见下图：

其中，(xpc, ypc)表示了计算域的参考原点的坐标（单位：m），alpc表示了计算网格与地理坐标系（笛卡尔坐标系）之间的偏转角；xlenc和ylenc分别表示了计算区域的长度和宽度（单位：m）；mxc和myc分别表示计算区域 x 和 y 方向上的网格数。对于本算例，我们可以看到计算域坐标原点为空间中的(0.0, 0.0)点，计算域坐标系和空间坐标系的方向重合；由于是水平向的一维流动，故ylenc和myc分别设定为 0.0 和 0；而计算域 x 方向的总长为 240.0 m，共包含360个网格，即 Δx = 2/3 m。
在垂直方向上，SWASH模型采用σ坐标系。VERT指令指定了垂向网格数，其必要参数为kmax，表示垂向总的网格数。若垂向各个σ层不等距，可输入thickness参数。thickness参数有两种输入形式，一是给定每一层的高度，单位为m，需要在指令尾加上M标识；另一种则给定每一个σ层所占总水深的百分比，可在指令尾加上PERC标识或省略（此为默认设置）。
需要注意的是，thickness参数是一连串的数字，各个数字间用空格隔开。各数字的顺序对应是从水面到水底。
示例：

VERT 20   10. 10. 10. 10. 7.5 7.5 6. 6. 5. 5. 4. 4. 3. 3. 2. 2. 1.5 1.5 1. 1.

后面的INPgrid BOTTOM 和 READinp BOTTOM 指令用于生成和读取水深（地形）数据。生成区域为一个矩形，这个过程类似于计算域网格的生成。其中 INPGRID BOTTOM 后面的参数如下所示：

与之前的xpc和ypc相似(xpinp, ypinp)表示读入地形数据矩形区域的原点，alpinp为该区域与地理笛卡尔坐标系的夹角。把这个矩形区在长宽方向上分成若干个等距的小网格，mxinp和myinp是矩形区两个方向的格子数字，dxinp和dyinp分别是x和y方向上小网格的大小（单位：m）。在本算例中，由于计算域在水平上是一维的，故所需的地形数据也可是一维的（一条线），故myinp和dyinp都为零。注：后面的可选参数在都没有在该模型中涉及，故不在本blog中解释。
READINP BOTTOM指令会让模型从文件中’a14prwav.bot’中读取地形数据。由于mxinp=1，则所需地形数据只需要两个，即需要这个网格最左端和最右端的底面深度。a14prwav.bot中的数据如下所示：

 1.0
 1.0

它表示计算域中，(0,0)处的静止水深为1.0m，(240m,0)处的静止水深也为1.0。其余计算域的网格点静止水深通过这两个点的数值进行线性插值，即全局的底面深度都为1.0m（即计算域为平底，底高程为 -1.0m）。

初始条件与边界条件

INIT zero
$
BOU SIDE W CCW BTYPE WEAK SMOO 1.7896 SEC CON REG 0.02 3.5791
BOU SIDE E CCW BTYPE RADIATION
SPON RI 40.

INIT zero 表示了初始的水位为0m，且初始的速度、紊动能和耗散率等均为0。
BOUnd 是设定边界条件的前缀，在本算例中，我们需要定义 x = 0 （WEST）和 x = 240m （EAST）两个边界。


首先对于West边界，这是一个波浪入射边界。后续的第一个参数CCW表示在west边界上，定义的方向是逆时针；但对于该算例，水平方向是是一维的，故这个参数并不会起作用。之后的BTYPE及后续参数表示了边界的类型；在本算例中BTYPE WEAK表示边界时弱反射类型，需要后续用傅里叶级数或时间序列来补充，一般用于波浪入射边界。之后的SMOO表示SMOOTHING，即使用此选项侯，将应用斜坡函数以平稳地载入边界条件；后续的 1.7896 SEC 表示这个斜坡函数的作用时间为1.7896秒。后面的CON REG 0.02 3.5791表示了输入的波浪参数是固定值（CONstant），波浪是规则波（REGular），且波高 h = 0.02m，波周期 per = 3.5791秒。
注：最后的[dir]表示入射方向，因为本算例是个垂向二维的流动，入射方向已经固定，无需再给定dir这个参数。

对于EAST边界，模型采用了RADIATION参数，即SOMMerfeld辐射边界。
注：在模型中，RADIATION和SOMMerfeld等同，这在手册中并未注明；但笔者推荐使用SOMMerfeld。

除了上述两个边界，模型还在计算域右侧设置了一个海绵层（SPONgelayer指令）。后面的参数RIght表示了边界层在计算域的右端，且厚度 width = 40.0m。

注：这里的参数RI和East等同；建议在自己建模时，统一使用East这种形式。

数值方法

NONHYDrostatic BOX PREC ILU
$
DISCRET UPW NONE
DISCRET UPW UMOM V NONE
DISCRET CORRDEP NONE
$
TIMEI 0.2 0.5

首先，模型是非静压模拟，所以要指定参数 NONHYDrostatic。其垂向的数值离散格式采用了Keller-box格式（该格式常用于波浪模型中）。此外，对于压力方程的求解，SWASH采用了ILU预处理算子。

DISCRETization为首的指令给定了动量方程的数值离散格式（主要是对流过程的离散格式）。以上各个指令的含义分别是：

1. UPW NONE：此指令省略了UMOM H，它表示了模型对于水平向（Horizontal）动量求解，**不采用（NONE）**迎风格式（UPWind），而采用标准中心差分格式；
2. UPW UMOM V NONE：表示了模型对于垂直向（Vertical）动量求解，**不采用（NONE）**迎风格式（UPWind），而采用标准中心差分格式；
3. CORRDEP NONE：对于速度计算点上的水深、水位等数据，模型中**不采用（NONE）**高阶的插值方法。
   
   
   在时间变化率的求解上，模型的参数通过TIMEI METH指令给定。默认采用显式格式（EXPLlicit）；后续参数cfllow和cflhig参数表示允许Courant数的最小和最大值。
   

输出设置

$
GROUP 'GAUGE' 90 90 1 1
TABLE 'GAUGE' NOHEAD 'a14prw01.tbl' TSEC WATL OUTPUT 000000.000 0.0125 SEC
$
TEST 1,0

GROUP指令定义了模型中的哪些点的位置将被输出。以上述指令为例，定义了名为“GAUGE”的一个网格点的组，这个组中包含了x坐标从90 ~ 90m且y坐标从1 ~ 1m的（也即x = 90m，y = 1m这一点）。注：由于计算域是垂向二维的，所以这里的iy1和iy2并不会真正地在模型中起作用。

之后我们将“GAUGE”中所有点的数据以数据表的格式输出到名为’a14prw01.tbl’的文件中。这个文件不包含文件头（NOHEADer），这样边于matlab的读取。这个文件中的时间数据是以秒为单位的（TSEC），输出的时间序列是水位（WATLEV）的时间序列。这个时间序列对应的起始时间点为00:00:00.000，时间间隔 deltblk = 0.0125秒（sec）。

执行参数

COMPUTE 000000.000 0.025 SEC 000112.000
STOP

COMPUTE指令可以让模型开始读入参数并运行。后面的参数依次表示了模型的起算时间点为 00:00:00.000，时间步为 0.025秒，模拟终止时间点为00:01:12.000。

模型运行

将五个参数文件a14prw01.sws、a14prw02.sws、a14prw03.sws、a14prw04.sws和a14prw05.sws，以及地形数据文件a14prwav.bot复制到swash.exe的同一目录下。并在这个目录下，用如下指令依次运行各个参数文件：

swashrun a14prw01

swashrun后面的带的是参数文件的文件名（不带扩展名）。
全部运行结束后我们会得到对应的5个.tbl文件和.prt文件。其中，.tbl文件即输出的结果数据；而.prt文件是模型运行过程中的输出。

结果数据的提取

数据可视化

将所有的.tbl文件复制到mkplot.m及mkplot01.m、mkplot02.m、mkplot03.m、mkplot04.m和mkplot05.m所在的文件夹中，并运行mkplot.m。上述的matlab脚本都在压缩包中。
mkplot.m实现了数据的可视化，于是得到以下的五个png文件。

1. prowave1.png
   
2. prowave2.png
   
3. prowave3.png
   
4. prowave4.png
   
5. prowave5.png
   

结果数据结构及matlab脚本原理

对于.tbl文件，其内容是两列数，第一列即时间（单位：秒），第二列是水位（单位：米）。在matlab中，我们可以通过load或textread指令进行读取；读取后，这些数据就能以一个两列的数组载入matlab中。
随后，我们需要对数据进行无量纲化，即所有的时间要表示为时间值与入射波周期的比值，所有水位要表示成水位值与入射波振幅的比值。
最后，利用plot指令将数据画成点线图，在利用 print -dpng 指令将点线图图窗导出为图片文件。

以mkplot01.m为例，matlab脚本如下：

a0=0.01;
T=3.5791;
tt=0:0.01:72;
w=2*pi/T;
y=cos(w*tt);
load a14prw01.tbl
t=a14prw01(:,1)/T;
wl=a14prw01(:,2)/a0;
plot(t(1:5:2881),wl(1:5:2881),'ko','MarkerSize',2);hold;
plot(tt/T,y,'k');
axis([10 20 -2 2])
xlabel('t/T_0');
ylabel('\zeta/a_0');
title('kh=\pi; 2 equidistant layers')
set(gca,'PlotBoxAspectRatio',[2 1 1]);
print -dpng prowave1.png

---

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

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