非静压模型NHWAVE学习（6）——波浪模拟算例学习（Periodic wave over a submerged bar）

本blog介绍了NHWAVE模型自带算例Periodic wave over a submerged bar的编译运行；并分享了该算例配置文件、模拟结果的学习记录。
关于该算例的详细描述见NHWAVE手册或Ma等人的论文 ¹。

此外，本人在NHWAVE原版源代码基础上，更新了垂向网格分层百分比可指定、被动示踪剂输移、温度场输入输出等模块，并修正了植被水流模型中的一些BUG。有兴趣的朋友请上我的GitHub下载代码。

欢迎各位朋友对此提出指导意见；大家可以一起交流，共同进步！

算例简介

该算例模拟的是一列周期波经过一个潜坝时的变化。这个算例也是广泛地被应用于模型性能的测试。
在NHWAVE中，我们建立一个垂向二维的数值水槽来模拟该算例。其底部形状及波浪入射的位置如下图所示。波浪从左侧入射，经过这个底部的坡后会发生浅化与变形。为了观察波浪形状与变形情况，我们设置了如下图所示的 a~f 共6个测点。

该模型的配置可在 /examples/submerged_bar 中找到。配置文件共有三个，分别是input.txt、depth.txt和stat.txt。

网格与地形

本算例的网格配置见input.txt的DIMENSION、GRID和VERTICAL GRID OPTION部分。

! --------------------DIMENSION---------------------------------
! cell numbers
Mglob = 1750
Nglob = 1
Kglob = 3
! ------------------------GRID----------------------------------
! grid sizes
DX = 0.02
DY = 0.02
! ---------------------VERTICAL GRID OPTION--------------------
! IVGRD = 1: uniform; 2: exponential
IVGRD = 1
GRD_R = 1.1

如上所示，计算区域水平向被平分成1760个网格，垂向被均分为3个网格。水平向网格的尺寸为0.02m。
模型的地形配置见input.txt的BATHYMETRY部分及depth.txt文件。

! ----------------------BATHYMETRY---------------------------
! if analytical bathymetry, set ANA_BATHY = T
! DEPTH_TYPE = CELL_CENTER if the water depth is defined at
! cell center, otherwise, DEPTH_TYPE = CELL_GRID 
DEPTH_TYPE = CELL_CENTER
ANA_BATHY  = F
DepConst = 0.3

该部分参数设置表示，读入地形点数据将被定义在网格中心，即depth.txt中数据个数必须与Mglob保持一致。此外，还要注意的是，上图所示的底部高程z(m)与我们在depth.txt载入的水深呈相反数关系。例如，对于计算域最左侧的第一个网格，其底高程为-0.4m，则在depth.txt中，第一个被载入的网格水深为0.4（默认单位是m）。

运行参数配置

本算例的核心是波浪模拟，故本blog重点学习波浪参数。在NHWAVE中，首先需要在边界上指定波浪的入射：

! -------------------BOUNDARY_TYPE--------------------------------
! bc_type=1: free-slip
!         2: no-slip
!         3: influx
!         4: outflux (specified eta)
!         5: bottom or wall friction
!         6: radiation bc
BC_X0 = 3
BC_Xn = 1
BC_Y0 = 1
BC_Yn = 1
BC_Z0 = 1
BC_Zn = 1

模型设定了左边界BC_X0为influx边界，即说明该边界上有波浪入射。其具体参数在WAVEMAKER部分中：

! ---------------------WAVEMAKER------------------------------
! wavemaker
! AMP - wave height; PER - wave period; DEP - incident water depth
! THETA - incident wave angle
! LEF_SOL - left boundary solitary wave, need AMP,DEP
! LEF_LIN - left boundary linear wave, need AMP,PER,DEP
! LEF_CON - left boundary cnoidal wave, need AMP,PER,DEP
! LEF_STK - left boundary stokes wave, need AMP,PER,DEP
! LEF_TID - left boundary tide wave, has to specify in subroutine
! LEF_JON - left boundary for JONSWAP spectrum
! RIG_LIN - right boundary linear wave, need AMP,PER,DEP,THETA
! INI_ETA - initial surface elevation specified in subroutine initial
! INT_LIN - internal wavemaker for linear wave
! INT_CON - internal wavemaker for cnoidal wave
! INT_SOL - internal wavemaker for solitary wave
! INT_JON - internal wavemaker for JONSWAP spectrum
! INT_SPC - internal wavemaker for 2D spectrum (need spc2d.txt)
! INT_IRR - internal wavemaker for irregular wave (need irr2d.txt)
! FLUX_LR - impose flux at both left and right boundaries
! FOCUSED - left boundary focusing wave packet (isolated whitecap)
! WAV_CUR - left boundary coexisting waves and currents
WAVEMAKER = LEF_LIN
AMP = 0.02
PER = 2.02
DEP = 0.40
THETA = 0.0
CUR = 0.0
sd_return = 0.0

模型的左边界采用线性波边界条件，其波高设置为0.02 m，波周期为2.02 s，水深为0.04 m。以上的桑格参数对应了上述代码中的 AMP，PER 和 DEP，它们设置线性波边界条件所必须的。

对于模型的监测点，我们需要在input.txt的PROBE OUTPUT部分中指定：

! --------------------PROBE OUTPUT---------------------------------
! output variables at stations which are given in file stat.txt
! in form of x,y,z, where z=-1.0 outputs velocity of all layers
NSTAT = 6
PLOT_INTV_STAT = 0.01

上述内容表示，一共有六个监测点，输出结果的间隔为0.01s；它们的位置详见stat.txt。

模型的编译运行及模拟结果

首先，读取结果文件夹中的 probe_[number] 文件（这里的number即监测点的编号）；其中第一列数据为时刻（s），第二列为自由水面高程（m）。
对比模型计算与物理实验中a~f点水位变化过程，如下图所示。图中的圆圈为物理实验结果，实线为模拟结果。可以看出，NHWAVE模拟结果与实验的拟合程度十分高。

---

1. Ma G , Shi F , Kirby J T . Shock-capturing non-hydrostatic model for fully dispersive surface wave processes[J]. Ocean Modelling, 2012, 43-44(22-35):22-35. ↩︎