SUNTANS模型学习（6）——学习lock-exchange算例

简介

Lock-Exchange是自然流体密度流或重力流的一个经典算例。在Lock-exchange流发展的过程中，不同密度流体的交界面会产生Kelvin–Helmholtz（K-H）不稳定，我们会观察到Kelvin–Helmholtz涡。
在静压模型中，我们无法模拟得到K-H不稳定性和K-H涡。故该算例可用于测试SUNTANS的非静压模拟能力。不过，代码包自带Lock-exchange算例的网格和参数并不能模拟出明显的K-H涡；想要模拟出明显的K-H不稳定性，详见Fringer等¹的论文。
对于Lock-Exchange实验，在初始时刻，我们将不同密度的两种水体分别放置在容器的左右两侧，并保持它们交界面的垂直。
该算例文件夹（/examples/lock-exchange），包含了两套计算域/网格配置：

1. 垂向二维的计算域；
2. 三维计算域。
   可在/rundata/suntans.dat中修改pslg参数来实现不同网格的使用。若pslg oned.dat，则采用垂向二维模拟；若pslg twod.dat，则采用三维模拟。

网格配置

上述两套网格配置采用了不同的平面网格布置，但采用了相同的垂向网格布置。
对于垂向二维的计算域，其平面网格布置为一串三角形，如下图所示：

对于三维的计算域，其平面网格布置如下图所示：

而两者的垂向采用上疏下密的网格布置（如下图所示），网格的拉伸比为1.025。垂向网格布置的参数在 /rundata/suntans.dat 中设定，即 rstretch = -1.025。

因此，两套网格的三维图分别如下：

1. 垂向二维网格
   
2. 三维网格
   

定解条件设置

初始条件设置

在本算例中，需要设定的初始条件为水深和初始盐度。相关设置都在 initialization.c 文件中。
初始水深通过函数 ReturnDepth 设定：

REAL ReturnDepth(REAL x, REAL y) {
  REAL length, xmid, shelfdepth, depth;
  return 20;
}

即初始水深为20m。此外，根据水深和垂向网格的拉伸比 rstretch，函数 GetDZ 确定垂向网格的布置：

int GetDZ(REAL *dz, REAL depth, REAL localdepth, int Nkmax, int myproc) {
  int k, status;
  REAL z=0, dz0, r = GetValue(DATAFILE,"rstretch",&status);

  if(dz!=NULL) {
    if(r==1) 
      for(k=0;k<Nkmax;k++)
	dz[k]=depth/Nkmax;
    else if(r>1 && r<=1.1) {    
      dz[0] = depth*(r-1)/(pow(r,Nkmax)-1);
      if(VERBOSE>2) printf("Minimum vertical grid spacing is %.2f\n",dz[0]);
      for(k=1;k<Nkmax;k++) 
	dz[k]=r*dz[k-1];
    } else if(r>-1.1 && r<-1) {    
      r=fabs(r);
      dz[Nkmax-1] = depth*(r-1)/(pow(r,Nkmax)-1);
      if(VERBOSE>2) printf("Minimum vertical grid spacing is %.2f\n",dz[Nkmax-1]);
      for(k=Nkmax-2;k>=0;k--) 
	dz[k]=r*dz[k+1];
    } else {
      printf("Error in GetDZ when trying to create vertical grid:\n");
      printf("Absolute value of stretching parameter rstretch must  be in the range (1,1.1).\n");
      exit(1);
    }
  } else {
    r=fabs(r);
    if(r!=1)
      dz0 = depth*(r-1)/(pow(r,Nkmax)-1);
    else
      dz0 = depth/Nkmax;
    z = dz0;
    for(k=1;k<Nkmax;k++) {
      dz0*=r;
      z+=dz0;
      if(z>=localdepth) {
	return k;
      }
    }
  }
}

上面的 if 结构表明，rstretch的大小只能在1~1.1之间，其正负号代表了垂向网格加密的位置，即正数表示顶部网格加密，负数表示底部网格加密。
初始盐度场通过函数 ReturnSalinity 设定：

REAL ReturnSalinity(REAL x, REAL y, REAL z) {
  return -.001*tanh(2.0*2.6467/5*(x-50.0));
}

上述表达式即：
$${\frac{\rho }{\rho }}_0=-\frac{1}{2}\frac{\Delta \rho }{{\rho }_0}tanh[\frac{2tan{h}^{-1}(\alpha )}{\delta }(x-\frac{L}{2})]$$
式中，Δρ/ρ₀=2×10^-3表示密度差，α=0.99，δ=5m 表示密度界面厚度，L=100m 表示计算域宽度。初始密度场如下图所示：

边界条件设置

在本算例中，除自由水面外，其余边界均为type=1（closed boundaries），即垂向速度为零的壁面边界。此时，边界条件无需通过 boundaries.c 来指定。

密度状态方程

在本算例中，水体密度仅由盐度s确定，其状态方程表达式体现在 state.c 文件的 StateEquation 函数中：

REAL StateEquation(const propT *prop, const REAL s, const REAL T, const REAL p) {
  return prop->beta*s;
}

其中，beta为定值参数，在 /rundata/suntans.dat 中指定。

其它参数配置

该算例的时间步长为Δt=0.2s （suntans.dat: dt = 0.2），总共运行500个时间步（suntans.dat: nstep = 500）。
考虑到网格尺度Δx = L/200，Voronoi点到三角形边的距离为D_g=Δx/sqrt(3) = 0.2887。水平速度的最大值约等于 u_max=0.5m/s（不会超过0.44m/s），垂向速度的最大值约为w_max=0.2m/s。所以Courant数C_x=u_maxΔt/D_g=0.35、C_x=w_maxΔt/Δz_min=0.2*0.2/0.78=0.05。
动量平流项采用中心差分格式（suntans.dat: nonlinear = 2），所以水平网格的Peclet数满足 Pe < 2/C_x。因此，对于水平稳定性，这要求ν_H ≥u²_maxΔt/2=0.025 m²s⁻¹ (suntans.dat: nuH = 0.025)。同样，对于垂直稳定性，ν_V= 0.016 m²s⁻¹ (suntans.dat: nu = 0.016)。

模拟结果

以下展示模拟得到的密度分布和速度矢量图。

参考文献

---

1. Fringer O B , Gerritsen M , Street R L . An unstructured-grid, finite-volume, nonhydrostatic, parallel coastal ocean simulator[J]. Ocean Modelling, 2006, 14(3):139-173. ](https://doi.org/10.1016/j.ocemod.2006.03.006) ↩︎