非静压模型SWASH学习（7）——自制算例Lock-Exchange

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)²
3. An accurate and efficient finite-difference algorithm for non-hydrostatic free-surface flow with application to wave propagation³

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

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

算例简介

Lock-Exchange 问题指的是一个封闭水槽中，由垂直隔板隔开的不同密度静止流体在去除隔板后的流动与混合问题。在去除隔板后，两个不同密度的流体在斜压作用下，分别在流体表层与底层形成两股等速但反向的重力流。在此过程中，不同密度流体的交界面会因剪切失稳而产生涡旋，影响了界面处的水体混合。
总之， Lock-Exchange 问题能描述重力流的诸多关键物理过程，是模拟重力流时空变化的经典算例。本博客将以 Fringer & Street (2006)⁴ 和 Lai et. al (2010)⁵ 的非静压模型研究为参考，尝试Lock-Exchange的静压模拟和非静压模拟，以展现两者的显著区别。

参考文献：

1. An unstructured-grid, finite-volume, nonhydrostatic, parallel coastal ocean simulator (Fringer & Street).
2. A nonhydrostatic version of FVCOM: 1. Validation experiments (Lai et. al).

模型输入文件及相关数据可下载于【传送门】。若要运行此算例，请将此下载的压缩文件解压，将含有.sws、 .bot 和 .sed 后缀的所有文件解压至swash.exe所在的目录下。

模型配置

模型的初始状态如下图所示。其中，左侧为低密度流体，右侧为高密度流体。在本文中，我们定义密度仅和水体盐度相关；左侧水体为清水，右侧水体为含盐量5.0ppt的盐水。

注：上图来自文献 A nonhydrostatic version of FVCOM: 1. Validation experiments (Lai et. al)。

网格及参数设置

网格与地形

MODE DYN ONED
$
CGRID 0. 0. 0. 0.8 0. 400 0
$
VERT 100
$
INPGRID BOTTOM 0. 0. 0. 1 0 0.8 0.
READINP BOTTOM 1. 'lock.bot' 1 0 FREE

由于模型采用了立面二维的计算区域。故在水平面上，模型采用一维模式（MODE DYN ONED）。将长0.8 m的计算域均分为400个网格，每个网格的长度 Δx = 0.002 m。此外，在垂直方向上，计算域被均分为100个σ网格。
随后，通过 INPGRID 和 READINP 读取了 lock.bot 文件中初始水深（即地形）数据。

初始条件与边界条件

INIT ZERO
$
INPTRANS  SEDIMENT 0. 0. 0. 400 0 0.002 0.
READTRANS SEDIMENT 1. 'lock.sed' 1 0 FREE

模型的初始水位、流速值均设为零（INIT ZERO）。而对于盐度，本模型通过READTRANS指令将 lock.sal 文件中的数据读入。
lock.sal 文件包含了401个数据，它们对应了401个格点上的盐度值。前一半的浓度值为 0.0，表示了低密度流体中不含泥沙；后一半的值为 5.0。

物理参数设置

VISC 0.000001
TRANSPORT 0.000001

设定模型的水平粘度（分子粘度）为10^-6m²/s，盐度的分子扩散系数为10^-6m²/s。

数值求解方法

NONHYDrostatic standard
$
$ Note that energy is conserved
DISCRET UPW HEAD
DISCRET UPW UMOM V NONE
DISCRET UPW WMOM H NONE
DISCRET UPW WMOM V NONE
DISCRET TRANSP H NONE
DISCRET TRANSP V NONE

与之前的波浪模拟不同，此处还要对物质输运方程的数值方法进行设定。设定指令为 DISCRET TRANSP。在本模拟中，对于盐度水平对流和垂直对流，采用中心差分格式进行离散，而非采用迎风格式（NONE）。

此外，为了表现非静压模拟与静压模拟的差异，我制作了两套参数文件，分别是 lock.sws 和 locksta.sws。两者间的主要区别就是 lock.sws 包含了 NONHYDrostatic 这一指令，而 locksta.sws 中不包含这个指令。当不包含 NONHYDrostatic 指令时，模型进行静压模拟计算。

模型输出

BLOCK 'COMPGRID' NOHEAD 'lock01.mat' LAY 3 XP ZK TSEC SEDK VELK OUTPUT 000000.000 1 SEC

模型的输出文件采用 .mat 格式，输出数据包含各网格x坐标值、各σ层的厚度、时间、泥沙浓度场和流速场。第一个输出数据对应了 000000.000 时刻，之后每隔1.0s输出一组数据。

**注：**对于参数文件 locksta.sws ，输出文件的名称为 lock02.mat。

计算时间

COMPUTE 000000.000 0.001 SEC 000030.000

设置模拟时间为30s；计算的时间步长 Δt = 0.0005 s。

模拟结果

将参数文件 lock.sws、locksta.sws，以及数据文件 lock.bot、lock.sal 复制到swash.exe的同一目录下。并在这个目录下，用swashrun分别运行 lock 和 locksta。

运行完成后，得到 lock01.mat 和 lock02.mat。
通过脚本mkplot_LE.m 将结果数据绘制成云图。matlab脚本文件也在压缩包中。得到云图如下所示。

1. t = 5.0 s
   
2. t = 10.0 s
   
3. t = 15.0 s
   
   上述结果展现了在重力流模拟中，非静压模拟和静压模拟的显著差异。静压模拟无法捕捉不同密度流体间的界面涡结构。

---

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.1002/fld.595 ↩︎

4. https://doi.org/10.1016/j.ocemod.2006.03.006 ↩︎

5. https://doi.org/10.1029/2009JC005525 ↩︎