非静压模型NHWAVE学习（14）—— 算例制作：开闸式异重流（lock-exchange flow）

NHWAVE学习—— 算例制作：开闸式异重流（lock-exchange flow）

算例简介
模型配置
代码修改及输入文件制作
- 代码修改
- 参数文件制作（input.txt）
- 水深和初始密度场文件制作（depth.txt & sali0.txt）
模型的编译运行及模拟结果
- 模型编译
- 结果可视化

本blog介绍了基于NHWAVE模型制作一个经典算例 lock-exchange flow，及其的编译运行；并分享了该算例配置文件、模拟结果的学习记录。
关于该算例的详细描述见相关论文¹^, ²。

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

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

最后，模型的相关代码及输入文件可下载于【传送门】。

算例简介

异重流又称为重力流、密度流，是自然界中常见的现象，例如河口的盐水入侵、深水湖泊中的潜流、含沙水流等。最典型的异重流为开闸式异重流（lock-exchange flow），即指封闭水槽中由隔板隔开的不同密度静止液体在去除隔板后的混合问题。

在去除隔板后，由于密度不同导致的斜压力的作用，在液体表层和底层形成了速度相同、方向相反的两股重力流，此过程比较复杂，涉及一些独特的物理现象:如盐度界面的盐淡水混合、Kelvin-Helmholtz 不稳定性导致的涡的形成等。正因如此，LE问题一直作为研究重力流时空变化的经典算例，从物理实验、数值模拟和理论研究方面被很多研究者所关注。³

该算例的实验装置和初始状态如下图所示。水槽被两种不同密度的流体充满，中间用隔板隔开。左半为低密度淡水，右半为高密度盐水（盐度为1.3592psu）。
在这里插入图片描述
模型中流体密度 $\rho$ 和盐度 $s$ 满足如下关系：
$\rho=999.972 \times (1+0.75 \times 10^{-3} s)$

模型配置

模型域为一个长0.8m，初始水深为0.1m的垂向二维平底水槽。计算区域的底面、侧壁均设置为可滑移的壁面边界（free-slip）。我们设置水平方向的空间步长为 $\Delta x= 2.0\times10^{-3}m$ ，网格数为400；垂向设置均匀的100个 $\sigma$ 层。

初始盐度分布满足：
$\begin{cases} 1.3592 &, x<0.4m,\\ 0.0 &, x > 0.4m \end{cases}$

此外，对于本次模拟，考虑无粘流（ $\nu=0$ ）和固定粘度流（ $\nu=1.0\times10^{-6} m^2/s$ ）两种情况。为了保持稳定，模拟过程中的CFL最大值设置为0.5。

代码修改及输入文件制作

代码修改

由于本算例中的密度-盐度方程并未内置于原来的代码中，所以我们需要将密度方程植入源代码的nhwave.F文件中。
（注：本次修改的基底代码为Kirby等人开发的代码，可从GitHub – JimKirby / NHWAVE 网站下载）

首先，打开nhwave.F文件，找到计算水体密度的子过程subroutine eval_dens。我们重点关注下面第30至42行代码，这部分代码即通过盐度（SALINITY）来计算密度Rho（这是一个数组）。

    subroutine eval_dens
!---------------------------------------------------------------------
!
!   equation of state
!
!---------------------------------------------------------------------use globalimplicit noneinteger, parameter :: kbb = 101integer, parameter :: kbbm1 = kbb-1real(SP), dimension(kbbm1) :: phy_z,rhoztmp,rhomeanreal(SP), dimension(Mloc,Nloc,kbbm1) :: rhozreal(SP), dimension(Kloc) :: zm,rhosinteger,dimension(1) :: reqreal(SP),dimension(:,:),allocatable :: xx,rhozlocreal(SP),dimension(Mglob,Nglob,kbbm1) :: rhozglob
# if defined (PARALLEL)integer,dimension(MPI_STATUS_SIZE,1) :: statusinteger,dimension(NumP) :: npxs,npys
# endifreal(SP) :: DELTZ,TF,SF,RHOF,HMAX,ETAMAXreal(SP) :: tmp1,tmp2,myvarinteger :: I,J,K,isum,jk,iglob,jglob,kk,n,len,nreq,NKloc! calculate density from equation of stateDO K = 1,KLOCDO J = 1,NLOCDO I = 1,MLOCIF(Mask(I,J)==0) cycle
# if defined (SALINITY)TF = Temp(I,J,K)SF = Sali(I,J,K)RHOF = SF*SF*SF*6.76786136E-6_SP-SF*SF*4.8249614E-4_SP+ &SF*8.14876577E-1_SP-0.22584586E0_SPRHOF = RHOF*(TF*TF*TF*1.667E-8_SP-TF*TF*8.164E-7_SP+ &TF*1.803E-5_SP)RHOF = RHOF+1.-TF*TF*TF*1.0843E-6_SP+TF*TF*9.8185E-5_SP-TF*4.786E-3_SPRHOF = RHOF*(SF*SF*SF*6.76786136E-6_SP-SF*SF*4.8249614E-4_SP+  &SF*8.14876577E-1_SP+3.895414E-2_SP)RHOF = RHOF-(TF-3.98_SP)**2*(TF+283.0_SP)/(503.57_SP*(TF+67.26_SP))Rho(I,J,K) = Rho0+RHOF
# endif# if defined (SEDIMENT)Rho(I,J,K) = (1.0-Conc(i,j,k))*Rho0+Conc(i,j,k)*SRho
# endifENDDOENDDOENDDO! find maximum water depthtmp1 = -largetmp2 = -largedo j = 1,Nlocdo i = 1,Mlocif(Mask(i,j)==0) cycle if(hc(i,j)>tmp1) tmp1 = Hc(i,j)  if(eta(i,j)>tmp2) tmp2 = Eta(i,j)enddoenddo
# if defined (PARALLEL)call MPI_ALLREDUCE(tmp1,myvar,1,MPI_SP,MPI_MAX,MPI_COMM_WORLD,ier)hmax = myvarcall MPI_ALLREDUCE(tmp2,myvar,1,MPI_SP,MPI_MAX,MPI_COMM_WORLD,ier)etamax = myvar
# endif! interpolate into physical z levelsdeltz = (hmax+etamax)/float(kbbm1)do k = 1,kbbm1phy_z(k) = (float(k)-0.5)*deltz-hmaxenddorhoz = Rho0do i = 1,Mlocdo j = 1,Nlocif(Mask(i,j)==0) cycledo k = 1,Kloczm(k) = sigc(k)*D(i,j)-Hc(i,j)    rhos(k) = Rho(i,j,k)enddocall sinter(zm,rhos,phy_z,rhoztmp,Kloc,kbbm1)do k = 1,kbbm1rhoz(i,j,k) = rhoztmp(k)enddoenddoenddo# if defined (PARALLEL)call MPI_GATHER(npx,1,MPI_INTEGER,npxs,1,MPI_INTEGER,  &0,MPI_COMM_WORLD,ier)call MPI_GATHER(npy,1,MPI_INTEGER,npys,1,MPI_INTEGER,  &0,MPI_COMM_WORLD,ier)NKloc = Nloc*kbbm1! put the data in master processor into the global var                                                           if(myid==0) thendo k = 1,kbbm1do j = Jbeg,Jenddo i = Ibeg,Iendiglob = i-Nghostjglob = j-Nghostrhozglob(iglob,jglob,k) = rhoz(i,j,k)enddoenddoenddoendifallocate(rhozloc(Mloc,NKloc))allocate(xx(Mloc,NKloc))do k = 1,kbbm1do j = 1,Nlocdo i = 1,Mlocjk = (k-1)*Nloc+jrhozloc(i,jk) = rhoz(i,j,k)enddoenddoenddo! collect data from other processors into the master processor                                                   len = Mloc*NKlocdo n = 1,NumP-1if(myid==0) thencall MPI_IRECV(xx,len,MPI_SP,n,0,MPI_COMM_WORLD,req(1),ier)call MPI_WAITALL(1,req,status,ier)do k = 1,kbbm1do j = Jbeg,Jenddo i = Ibeg,Iendiglob = npxs(n+1)*(Iend-Ibeg+1)+i-Nghostjglob = npys(n+1)*(Jend-Jbeg+1)+j-Nghostjk = (k-1)*Nloc+jrhozglob(iglob,jglob,k) = xx(i,jk)enddoenddoenddoendifif(myid==n) thencall MPI_SEND(rhozloc,len,MPI_SP,0,0,MPI_COMM_WORLD,ier)endifenddodeallocate(rhozloc)deallocate(xx)if(myid==0) thenrhomean = zerodo k = 1,kbbm1isum = 0do j = 1,Nglobdo i = 1,Mglobif(-HCG(i,j)<=phy_z(k)) thenisum = isum+1rhomean(k) = rhomean(k)+rhozglob(i,j,k)endifenddoenddoif(isum>=1) thenrhomean(k) = rhomean(k)/float(isum)elserhomean(k) = rhomean(k-1)endifenddoendifcall MPI_BCAST(rhomean,kbbm1,MPI_SP,0,MPI_COMM_WORLD,ier)# elserhomean = zerodo k = 1,kbbm1isum = 0do j = Jbeg,Jenddo i = Ibeg,Iendif(-Hc(i,j)<=phy_z(k)) thenisum = isum+1rhomean(k) = rhomean(k)+rhoz(i,j,k)endifenddoenddoif(isum>=1) thenrhomean(k) = rhomean(k)/float(isum)elserhomean(k) = rhomean(k-1)endifenddo# endif! linearly interpolate to obtain density at signa levelsRmean = Rho0do i = 1,Mlocdo j = 1,Nlocif(Mask(i,j)==0) cycledo k = 1,Kloczm(k) = sigc(k)*D(i,j)-Hc(i,j)enddocall sinter(phy_z,rhomean,zm,rhos,kbbm1,Kloc)          Rmean(i,j,1:Kloc) = rhosenddoenddoreturn  end subroutine eval_dens

为此，我们在第41行代码之后加上与我们密度方程对应的代码，以覆盖第33-41行的计算过程。修改后，这部分代码如下所示：

# if defined (SALINITY)TF = Temp(I,J,K)SF = Sali(I,J,K)RHOF = SF*SF*SF*6.76786136E-6_SP-SF*SF*4.8249614E-4_SP+ &SF*8.14876577E-1_SP-0.22584586E0_SPRHOF = RHOF*(TF*TF*TF*1.667E-8_SP-TF*TF*8.164E-7_SP+ &TF*1.803E-5_SP)RHOF = RHOF+1.-TF*TF*TF*1.0843E-6_SP+TF*TF*9.8185E-5_SP-TF*4.786E-3_SPRHOF = RHOF*(SF*SF*SF*6.76786136E-6_SP-SF*SF*4.8249614E-4_SP+  &SF*8.14876577E-1_SP+3.895414E-2_SP)RHOF = RHOF-(TF-3.98_SP)**2*(TF+283.0_SP)/(503.57_SP*(TF+67.26_SP))Rho(I,J,K) = Rho0+RHOF! for lock-exchange flowRho(I,J,K) = 999.972_SP * (1.0_SP + 0.75E-3_SP * SF )
# endif

所添加的代码中，SF表示所在点的盐度，这在原始代码的32行已定义了。

参数文件制作（input.txt）

网格设置
根据前面“模型配置”的描述，我们设置x方向的网格共400个，y方向的网格1个，z方向的 $\sigma$ 层为100（均匀分布）。

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

时间参数设置
本算例共运行30.0s，且每0.2s输出一次结果。每一次计算的时间步长通过CFL数控制，CFL数设置为0.5。

! -----------------------TIME----------------------------------
! time: total computational time/ plot time / screen interval 
! all in seconds
SIM_STEPS = 1000000000
TOTAL_TIME = 30.0
PLOT_START = 0.0
PLOT_INTV =  0.2
SCREEN_INTV = 0.05
! --------------------COURANT_NUMBER---------------------------------
CFL = 0.5

初始条件设置
初始水深（地形）通过depth.txt读入，它的数据定义在网格中心，与其一一对应；因此，需设置DEPTH_TYPE = CELL_CENTER。对于初始盐度场，我们通过sali0.txt文件读入；而初始温度场不需要设置，故INITIAL_EUVW = F、INITIAL_SALI = T。

! ----------------------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
! -------------------INITIAL CONDITION---------------------------
! if INITIAL_SALI = T, need file sali0.txt
INITIAL_EUVW = F
INITIAL_SALI = T

边界条件设置
根据前面“模型配置”中的描述，模型的壁面均为free-slip边界。所以，我们需要在input.txt文件中这样设置：

! -------------------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 = 1
BC_Xn = 1
BC_Y0 = 1
BC_Yn = 1
BC_Z0 = 1
BC_Zn = 1
! ----------------------BOTTOM ROUGHNESS-------------------
! Ibot = 1: given the drag coefficient Cd0
! Ibot = 2: given the bottom roughness height Zob
Ibot = 1
Cd0 = 0.000
Zob = 0.000
Dfric_Min = 0.0

物理参数设置
模型考虑了斜压效应，所以我们需要将正压模拟参数BAROTROPIC关闭，即BAROTROPIC = F。且模型考虑了动水压力，所以我们需要将非静压参数NON_HYDRO打开。

! ---------------------BAROTROPIC--------------------------
! if barotropic run, set BAROTROPIC = T
BAROTROPIC = F
! ----------------------NON-HYDRO---------------------------
! if non-hydrostatic simulation
NON_HYDRO = T

对于无粘流动模拟，我们设置 VISCOUS_FLOW = F 即可。而对于固定粘度流动模拟，我们需要如下设置：

! --------------------VISCOSITY------------------------------
VISCOUS_FLOW = T
IVTURB = 1
IHTURB = 1
PRODTYPE = 3
VISCOSITY = 1.e-6
Schmidt = 1.0
Chs = 0.0
Cvs = 0.0
RNG = T

其中IVTURB = IHTURB = 1表示垂向和水平向粘度均为固定值，且这个固定值分别为 VISCOSITY+Cvs 和 VISCOSITY+Chs。考虑到本算例的配置，其垂向与水平向粘度均为定值且相等，故我们设置 VISCOSITY = 1.e-6，Chs = Cvs = 0.0。

数值计算参数设置
首先，我们需要设置空间差分精度和时间差分的精度均为二阶，即HIGH_ORDER 和 TIME_ORDER 均为 SECOND。要注意的是，为了提高模拟精度，我们将启用 HLLC 求解器，它比默认的 HLL 求解器精度更高，更适合密度流的模拟。此外，还有一个需要注意的参数是干-湿临界水深，由于本算例不涉及干网格，故我们需要设置一个小于水深的MinDep值。

! ----------------------NUMERICS----------------------------
! Scalar convection scheme: "TVD" or "HLPA"
HIGH_ORDER = SECOND
TIME_ORDER = SECOND
CONVECTION = HLPA
HLLC = T
! ----------------------WET-DRY-------------------------------
! minimum depth for wetting-drying
MinDep = 0.005

输出结果设置
在本算例中，我们主要关心水深，x、z方向的流速以及盐度分布。对此我们需要将OUTPUT部分中的OUT_H、OUT_E、OUT_U、OUT_W、OUT_I设置为T。其余选项可以均设置为F。

此外，对于本算例中未涉及的模块，如植被模块、滑坡体模块等，我们需检查其是否关闭。例如，对于造波模块WAVEMAKER，我们最好将WAVEMAKER参数设置为 “WAVEMAKER = NONE”；这样，这个模块将不影响我们的计算。其余的设置详见我上传的压缩包中的input.txt文件。

水深和初始密度场文件制作（depth.txt & sali0.txt）

初始水深（depth.txt）的数据格式是一个共Nglob行，每行含Mglob个数字的数组。对于每一行，后面的数对应的网格在其前面数字的右侧（东侧）；靠下行的数据所对应的网格都在其上方行数据所对应网格的上侧（北侧）。

在本算例中，y方向仅一个网格，即Nglob=1；因此，depth.txt中应该只有一行数据，且这行数据共Mglob=400个。根据初始条件，这400个数字均为0.10。总结一下，depth.txt中包含一行数据，共400个0.10。

对于初始盐度场（sali0.txt），它包含了Kglob个有Nglob行、每行含Mglob个数的数组；且对应底部 $\sigma$ 层的数据在前，对应靠上层的数据在后。每个Nglob行、Mglob列的数组的数据格式与depth.txt中类似。

在本算例中，sali0.txt将包含Kglob=100行数，每一行包含Mglob=400个数据。根据初始条件，我们将每一行前200个数设置为0.0，后200个数据设置为1.3592。

对于水深和初始密度场文件制作，我在上传的压缩包中附了一个matlab脚本（InputFileMaker.m），它可以自动生成depth.txt和sali0.txt文件。

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

模型编译

在模型编译前，需要设置Makefile中的FLAG。对于该算例，我们需要启用SALINITY和TEMPERATURE模块（虽然TEMPERATURE模块与本算例无关，但是只要计算密度，不启用TEMPERATURE就会出错，这也是Kirby等人原始代码的一个问题）。

         FLAG_1 = -DDOUBLE_PRECISIONFLAG_2 = -DPARALLEL
#         FLAG_3 = -DLANDSLIDEFLAG_4 = -DSALINITYFLAG_5 = -DTEMPERATURE
#         FLAG_6 = -DBUBBLE
#         FLAG_7 = -DSEDIMENT
#         FLAG_8 = -DVEGETATION
#         FLAG_9 = -DINTEL
#         FLAG_10 = -DBALANCE2D
#	      FLAG_11 = -DOBSTACLE
#	      FLAG_12 = -DTWOLAYERSLIDE
#	      FLAG_13 = -DCORALREEF
#	      FLAG_14 = -DPOROUSMEDIA
#         FLAG_15 = -DFROUDE_CAP
#         FLAG_16 = -DCOUPLING
#         FLAG_17 = -DFLUIDSLIDE
#         FLAG_18 = -DLANDSLIDE_COMPREHENSIVE
#         FLAG_19 = -DDEFORMABLESLIDE

随后，通过make指令编译得到nhwave程序。之后，将input.txt、depth.txt和sali0.txt文件置于nhwave相同的目录下，并创建结果一个存放结果的文件夹results。最后，通过mpirun指令并行运行该程序。

结果可视化

在上传的压缩包中附了一个matlab脚本（ReadResults.m），它可以读取想要的时间步的结果文件，并画出密度分布云图和速度矢量图。

以无粘性流动（VISCOUS_FLOW = F）为例，下面即ReadResults.m生成的T = 5.0s、10.0s、15.0s和20.0s时刻的结果图。

T = 5.0s
T = 10.0s
T = 15.0s
T = 20.0s

Hu L, Xu J, Wang L, Zhu H. Effects of Different Slope Limiters on Stratified Shear Flow Simulation in a Non-hydrostatic Model. Journal of Marine Science and Engineering. 2022; 10(4):489. https://doi.org/10.3390/jmse10040489. ↩︎
Shin, J.O., Dalziel, S.B., Linden, P.F. Gravity currents produced by lock exchange. J. Fluid Mech. 2004, 521, 1–34. ↩︎
时健. 河口盐淡水垂向混合的非静压模拟研究. 博士论文. 2016. 南京. ↩︎